Bayesian model and selection signature analyses reveal risk factors for canine atopic dermatitis

Abstract

Canine atopic dermatitis is an inflammatory skin disease with clinical similarities to human atopic dermatitis. Several dog breeds are at increased risk for developing this disease but previous genetic associations are poorly defined. To identify additional genetic risk factors for canine atopic dermatitis, we here apply a Bayesian mixture model adapted for mapping complex traits and a cross-population extended haplotype test to search for disease-associated loci and selective sweeps in four dog breeds at risk for atopic dermatitis. We define 15 associated loci and eight candidate regions under selection by comparing cases with controls. One associated locus is syntenic to the major genetic risk locus (Filaggrin locus) in human atopic dermatitis. One selection signal in common type Labrador retriever cases positions across the TBC1D1 gene (body weight) and one signal of selection in working type German shepherd controls overlaps the LRP1B gene (brain), near the KYNU gene (psoriasis). In conclusion, we identify candidate genes, including genes belonging to the same biological pathways across multiple loci, with potential relevance to the pathogenesis of canine atopic dermatitis. The results show genetic similarities between dog and human atopic dermatitis, and future across-species genetic comparisons are hereby further motivated.

Fig. 1. BayesR and XP-EHH regions associated with canine AD in four dog breeds.

The results from BayesR analyses are presented for LR (a), GR (b), GSD (e), and WHWT (f). The red line defines the cutoff for AD effect variants at absolute effect size of 0.0001 resulting in 11 AD-associated loci in LR, one locus in GR, one in GSD, and two in WHWT above this threshold. The closest protein-coding gene(s) to the top effect variants are specified for each locus. Panels c, d and g, h present genome-wide scans for signatures of selection in canine AD cases or controls for each breed and candidate regions under selection are defined by selection variants above the threshold line (orange; -log₁₀(p) XP-EHH = 4). Protein coding gene(s) closest to the top selection variants for each candidate region of selection are highlighted in orange.

Fig. 2. Canine AD risk index in Labrador retrievers.

In LR, the risk genotypes (no risk alleles = 0, one risk allele = 0.5, and two risk alleles = 1 per locus) at the top effect variant for each of the 11 AD-associated loci were combined into a risk index. The distributions of risk indexes were shifted with higher values in cases compared to controls and the mean risk index was significantly higher in cases (mean = 15.6) compared to controls (mean = 13.1; two-sided t-test p = 1.52 × 10⁻²², t-statistic = 10.6, n = 321 dogs). The difference is also visualized with boxplots, red (cases) and blue (controls), indicating median, first and third quartiles, and range of whiskers defined by max 1.5 of box length (Supplementary Data 2).

Fig. 3. Canine AD-associated locus on chromosome 17 in Labrador retrievers.

One AD-associated locus in LR was located on chromosome 17 (a) with nine effect variants comprising a region of ~2 Mb. One of the top associated variants in human AD meta-GWAS (Sliz et al., 2021, purple) was located ~14 kb from effect variant chr17:c and one AD-associated SNP from the human AD multi-ancestry meta-GWAS (Paternoster et al., 2015, purple) was located in between variants chr17:h and chr17:i. Variants confirmed in sequenced dogs and in LD (r² > 0.8) with at least one of the effect variants extend >3 Mb (57.09–60.41 Mb), about 0.5 Mb from the epidermal differentiation complex (EDC) gene region (b). Black blocks show liver TADs in dogs. Black bars show 236 LD variants (plus two outside the figure boundaries at ~57.1–57.3 Mb) and 129 novel variants (plus four at ~57.3 Mb outside the figure boundaries). Variants with potentially regulatory functions are highlighted in light blue (overlap with canine ATAC-seq or both ATAC and cCRE), dark blue (ATAC and GeneHancer elements), and red (both ATAC, GeneHancer, and cCRE). Gray blocks show homozygous regions unique to the two ONT sequenced controls and not detected in the two sequenced cases (c). One homozygous region (56 kb) was found between effect variants chr17:b-c upstream of BCL9, and another region of homozygosity (14 kb) was spanning half of ACP6 (d). One 106 kb homozygous region, including six LD-variants and five novel variants, overlapped the VPS45 gene (e).

Fig. 4. A selection signal in Labrador retriever AD cases targets the TBC1D1 gene.

The top XP-EHH region in LR was located on chromosome 3. Information about Swedish LR kennels, from UK owner questionnaires, and information regarding Switzerland working dogs suggest that LR of gundog type are present in the low PC1 cluster whereas the common type LR are represented in the high PC1 cluster. This is further supported by the chocolate coat color only found among the common type LR (Supplementary Fig. 3). The low PC1 cluster is subsequently referred to the gundog subpopulation, whereas the common type LR are in the high PC1 cluster. A heavier body type is typically observed in the common type, whereas the working type LR are generally thinner as illustrated by photos of typical gundog and common LR types (a). A signal of selection was detected around 74.0-74.3 Mb on chromosome 3 overlapping with the TBC1D1 gene (associated with body weight). Two canine liver TADs (black bars) spanned the region. Red bars highlight variants with phyloP >2.56 (b-c). The allele C at the top variant position chr3:74,218,744 (chr3:sel) in LR was more frequent in the common type LR (a). A higher iHH in cases (618 kb) compared to controls (205 kb) was observed for the allele C at chr3:sel. The EHH plot of allele C indicates that LD decay is not complete until around 69 Mb and 77 Mb, respectively, in cases (d), and 72 Mb and 76 Mb in controls (e). Twelve genes were located within the candidate region under selection (black box). Association with canine AD within the extended region was calculated using plink --assoc (gray bars) and plink --logistic (covariates PC1 and PC2, black bars) resulting in associated variants along the haplotype; top variants (purple) named chr3:assocA through chr3:assocD (f). A high frequency (57.2%) was observed for the risk haplotype carrying risk alleles at chr3:assocA-chr3:assocD and allele C at chr3:sel, whereas nine different haplotypes (frequencies 1.0–9.9%) carried a mix of risk and non-risk alleles (g). LD was higher when considering only the chr3:assocA-chr3:assocC risk haplotype, CCG (frequency of 76.4%), compared to haplotype frequencies of 1.5–8.8% for the other combinations (h). Homozygous CCG was found in 129 (72.5%) of the cases and 68 (47.6%) of the controls (i; Supplementary Data 8).

Fig. 5. A selection signal in German shepherd controls spans the LRP1B gene.

Dogs from GSD kennels with a higher proportion of working compared to show merits and with gray or black coat color were more frequent in the low PC1 cluster whereas the opposite was seen for dogs in the high PC1 cluster (Supplementary Fig. 4). The subpopulations were subsequently referred to as the working (low PC1) and show (high PC1) type subpopulations. Dogs homozygous for the T allele were more frequent in the working type and the pictures illustrate the common gray color in a working type GSD and a show type with the typical GSD color pattern of brown with black saddle and a sloping cross (a). The selection signal on chr19 consisted of 1078 selection variants (-log₁₀(p) XP-EHH ≥4) spanning the LRP1B gene and with the top selection variant located at chr19:44,248,511 (chr19:sel). The KYNU gene was positioned within the same canine liver TAD as LRP1B (black blocks). Red bars indicate selection variants with high phyloP scores (>2.56; b, c). The iHH for chr19:sel for allele T in cases (3.18 Mb) was lower compared with controls (6.47 Mb) indicating selection in controls (d, e). The association with canine AD (plink --assoc model, gray bars) was high across LRP1B and overlapping with chr19:sel. The signal was lost in the logistic model (black bars) where PC1, PC2, -log₁₀(IgA), and -log₁₀(Age) were included as covariates. Protein coding genes are presented within the XP-EHH region (black box; f; genes outside this region and three small transcripts without official gene symbols (ENSCAFG00000024630, CFRNASEQ_PROT_00059550, and ENSCAFG00000030677) within the region were excluded from the plot).

Fig. 6. Allele frequency differences between breed types in both Labrador retrievers and German shepherds.

In LR, the number of cases with CCG (risk haplotype on chromosome 3) was significantly higher compared to controls of the common type (χ² = 17.5, p = 2.81 × 10⁻⁵) whereas there was no significant difference in the gundog type (χ² = 1.06, p = 0.303; homozygous = 2, heterozygous = 1, or homozygous non-risk = 0; a). In GSD, a homozygous T/T genotype at the top selection variant, chr19:sel, was more frequent among controls of working type (χ² = 5.21, p = 0.0224) but in the show type there was no difference (χ² = 0, p = 1; b). Dogs were split into subpopulations by a cutoff at PC1 of −0.05 in LR and 0 in GSD. Pearson’s Chi-squared test with Yates’ continuity correction of allele frequencies between cases and controls within each breed type was used. Two GSD controls from the working subpopulation did not have a genotype call and were excluded from the plot and this calculation (Supplementary Data 9).