A comprehensive analysis of gut and skin microbiota in canine atopic dermatitis in Shiba Inu dogs
Microbiome volume 11, Article number: 232 (2023)
Like its human counterpart, canine atopic dermatitis (cAD) is a chronic relapsing condition; thus, most cAD-affected dogs will require lifelong treatment to maintain an acceptable quality of life. A potential intervention is modulation of the composition of gut microbiota, and in fact, probiotic treatment has been proposed and tried in human atopic dermatitis (AD) patients. Since dogs are currently receiving intensive medical care, this will be the same option for dogs, while evidence of gut dysbiosis in cAD is still missing, although skin microbial profiling in cAD has been conducted in several studies. Therefore, we conducted a comprehensive analysis of both gut and skin microbiota in cAD in one specific cAD-predisposed breed, Shiba Inu. Additionally, we evaluated the impact of commonly used medical management on cAD (Janus kinase; JAK inhibitor, oclacitinib) on the gut and skin microbiota. Furthermore, we genotyped the Shiba Inu dogs according to the mitochondrial DNA haplogroup and assessed its association with the composition of the gut microbiota.
Staphylococcus was the most predominant bacterial genus observed in the skin; Escherichia/Shigella and Clostridium sensu stricto were highly abundant in the gut of cAD-affected dogs. In the gut microbiota, Fusobacteria and Megamonas were highly abundant in healthy dogs but significantly reduced in cAD-affected dogs. The abundance of these bacterial taxa was positively correlated with the effect of the treatment and state of the disease. Oclacitinib treatment on cAD-affected dogs shifted the composition of microbiota towards that in healthy dogs, and the latter brought it much closer to healthy microbiota, particularly in the gut. Additionally, even within the same dog breed, the mtDNA haplogroup varied, and there was an association between the mtDNA haplogroup and microbial composition in the gut and skin.
Dysbiosis of both the skin and the gut was observed in cAD in Shiba Inu dogs. Our findings provide a basis for the potential treatment of cAD by manipulating the gut microbiota as well as the skin microbiota.
Similar to humans, dogs suffer from spontaneously occurring diseases, such as cancers, and age-associated diseases [1,2,3]. Allergy is one such disease, and atopic dermatitis (AD) is the most common chronic inflammatory skin disorder both in humans and in dogs. In both species, AD is a complex and heterogeneous disease that is thought to be affected by multiple factors, such as host genetics and environmental factors [4,5,6]. In dogs, canine AD (cAD) exhibits similar clinical immunological characteristics to those of human AD . cAD has a strong association with breed, although the breed predisposition varies among geographical locations , supporting a complex interaction of gene and environmental factors. Environmental factors, such as pH, temperature, hydration, hygiene practices, and food, play a pivotal role in maintaining the homeostasis of the host-microbiome , and changes in microbial composition, namely, dysbiosis, are known to drive the pathology of complex diseases . In fact, dysbiosis of the skin microbiota has been the most well described in human AD patients, primarily defined as an overgrowth of Staphylococcus spp. on the human AD skin [9, 11]. Similarly, several studies have profiled skin microbiota in cAD-affected dogs, demonstrating skin dysbiosis in cAD [12,13,14]. Two of these studies have demonstrated a predominance of the genus Staphylococcus on the cAD dog skin compared to the healthy dog skin, confirming culture-based analysis of cAD skin that demonstrated the higher prevalence of this bacterial genus in cAD-affected dogs compared to healthy dogs .
Furthermore, a number of studies have demonstrated that gut dysbiosis is accompanied in human AD patients (reviewed in ). This is not surprising because both the gut and skin are abundantly colonised by distinct microbial communities, and these organs share a number of important anatomical and physiological characteristics . In fact, although their relevance is underestimated and the mechanism is often not fully understood, an association between gut and skin manifestations has been reported in other pathological conditions in patients (reviewed in ). To date, the role of the gut microbiota has been described as an immune system regulator by its production of bacterial metabolites. The most well-known example is short-chain fatty acids (SCFAs), including butyrate, propionate, and acetate, which have anti-inflammatory and immunomodulatory functions in other organs and the skin [18, 19]. As such, many studies suggest that probiotic supplementation may be a potentially effective treatment option for human AD; however, further studies are needed . This is considered to be a similar situation for the clinical management of cAD. While several studies of gut microbiota in cAD-affected dogs have recently been reported; a pilot study using a limited number of Beagles  and one breeding colony of Shiba Inu dogs , or a study of multiple dog breeds , no studies dissecting both skin gut microbiota in cAD-affected dogs are available to date.
This is the first study that describes a comprehensive microbiome profile of both skin and faeces obtained from cAD and healthy Shiba Inu dogs. This dog breed is highly susceptible to cAD in Japan, as shown by a high number of cases in published studies of cAD in this region [24,25,26]. More specifically, we analysed the microbiota composition of 12 skin sites, which were evaluated for clinical score (the fourth version of the Canine Atopic Dermatitis Extent and Severity Index; CADESI-04), as well as intestinal microbiota in faecal samples of healthy and cAD Shiba Inu dogs. In addition, we were able to evaluate the impact of Janus kinase (JAK) inhibitor oclasitinib (Apoquel®; Zoetis, Parsippany, NJ, USA) treatment, which is commonly used for cAD treatment, on the gut and skin microbiota of cAD-affected dogs. We focused on one specific dog breed to minimise the potential influence of the host nuclear genome on the skin and intestinal microbial community between samples. Meanwhile, the impact of the mitochondrial DNA (mtDNA) has not been evaluated in this dog breed, particularly in the context of microbiota as well as pathological relevance. Since our group and others have demonstrated that variations in mtDNA are associated with the composition of microbiota in humans and rodents [27, 28], we sequenced the whole mitochondrial genome of Shiba Inu dogs to evaluate the potential association between mtDNA variants and microbiota composition in cAD-affected dogs.
Sample subjects and sampling for microbiome analysis and mitochondrial genome sequencing
In total, 40 dogs were recruited, of which 20 were healthy dogs (Healthy). Shiba Inu dogs diagnosed with cAD at Animal Medical Centre, Faculty of Agriculture, Tokyo University of Agriculture and Technology and private veterinary practices in Japan were recruited for this study. All dogs with cAD fulfilled more than five of Favrot’s diagnostic criteria, and other potential causes for pruritis, such as ectoparasite infestations and bacterial or fungal infections were excluded, according to a detailed guideline . CADESI-04 scores for the same case were monitored by the same veterinary clinicians (TO and YS). The Pruritis Visual Analog Scale (PVAS) was graded by participating cAD-affected dog owners to assess the severity of pruritic manifestation . The clinical scores (CADESI-04 and PVAS) of cAD-affected dogs are presented in Figure S1.
Of the 20 cAD dogs, 10 dogs were newly diagnosed with cAD without any previous medical treatments (cADpre), and 10 were already diagnosed with cAD before the study cohort recruitment and had been on oclacitinib treatment (cADtreat; Figure S1a). The age and sex of the recruited dogs are summarised in Table S1. The ten dogs newly diagnosed with cAD (cADpre) were consequently under oclacitinib treatment for 2 weeks after the diagnosis. After the 2-week treatment, the dogs were clinically re-evaluated and skin swab and faecal samples were obtained (cADpost). cADtreat dogs had been on oclacitinib treatment for an average of 13.24 months (standard deviation, s.d. ± 8.458). Additionally, 20 healthy Shiba Inu dogs were recruited upon their visits for annual vaccination or general health checks. The dog groups tested in this study is depicted in Figure S1a. Clinical history, housing environment, and dietary information were collected and recorded (Supplementary Material S1).
Skin swab samples were obtained from the perilesional area of 12 body sites, which are included in the CADESI-04 evaluation. More specifically, these included perilabial area (lips), medial pinnae (concave pinnae), axilla, front paw (dorsal and palmar sides combined), hind paw (dorsal and palmar sides combined), cubital flexor (elbow fold), palmar metacarpal (from carpal to metacarpal pad), flank, inguinal area (groin), abdomen, perineum (from genitalia to anus) and ventral tail (proximal) (Figure S1b) . All dogs did not receive topical disinfectants or shampoos at least 3 days before the sampling. Samples were collected from the left side of the body sites unless skin lesions were present only on the right side; since there were no cases of unilateral right-sided lesions, all samples were collected from the left body side.
Sterile culture swab applicators were soaked in SCF-1 buffer (50 mM Tris, 1 mM EDTA, 0.5% Tween 20; Teknova, CA, USA) and rubbed on the perilesional skin 40 times, while the swab was rotated every 10 strokes. Swabs were stored in tubes filled with 400 µl of SCF-1 buffer at − 80 °C until bacterial DNA isolation was performed. Faecal samples were also collected from the dogs by inserting the swab into the rectum and stored in tubes filled with 400 µl SCF-1 buffer at − 80 °C. A total of 600 skin samples, 50 faecal samples and 40 buccal samples were collected from the study cohort. A summary of the sample strategy is displayed in Figure S1a.
Bacterial DNA isolation, library preparation and sequencing of the bacterial 16S rRNA gene
Skin microbial DNA isolation, library preparation and sequencing
Microbial DNA was extracted from the skin swab samples using a QIAmp DNA Microbiome Kit (Qiagen, Hildfeld, Germany) according to the manufacturer’s protocol. Skin microbial DNA samples were further processed for library preparation for the hypervariable V1–V2 region of the 16S rRNA gene, as previously described [28, 32, 33]. The final library was sequenced on the Illumina MiSeq platform using v2 chemistry, generating 2 × 250bp reads.
Faecal microbial DNA isolation, library preparation and sequencing
Microbial DNA was prepared from the faecal samples using a Qiagen Power Soil Kit (Qiagen) according to the manufacturer’s protocol. Faecal microbial DNA samples were further processed for library preparation for the 16S rRNA gene. The 16S rRNA gene was amplified using uniquely barcoded primers flanking the V3 and V4 hypervariable regions (341F - 806R) with fused MiSeq adapters and heterogeneity spacers in a 25 µl volume, consisting of one µl template DNA, four µl of each forward and reverse primer, 0.25 µl Phusion Hot Start II DNA polymerase (0.5 units), 0.5 µl dNTPs (200 µM each) and 5 µl of HF buffer. The PCR conditions were as follows: initial denaturation for 30 s at 98 °C; 32 cycles of 9 s at 98 °C, 30 s at 55 °C, and 45 s at 72 °C; and a final extension for 10 min at 72 °C. The concentration of the PCR products was estimated on 1.5% agarose gels using the image software Quantum Capt v16.04 (Vilber Lourmat Deutschland GmbH, Eberhardzell, Germany) with a known concentration of DNA marker (100 bp DNA ladder; Thermo Fisher Scientific GmbH, Dreieich, Germany) as the internal standard for band intensity measurement. The PCR products were pooled into approximately equimolar subpools, as indicated by band intensity, followed by clean-up of the subpools using a GeneJet column purification kit (Thermo Fisher). The concentrations of subpools were quantified using a Qubit dsDNA BR Assay Kit on a Qubit fluorometer (Thermo Fisher Scientific GmbH). The subpools were combined in one equimolar final pool for each library. The final pools were cleaned up with magnetic beads (MagSi-NGSPREP Plus; Steinbrenner Laborsysteme GmbH, Wiesenbach, Germany), the concentration was measured by qPCR using a NEBNext Library Quant Kit for Illumina (New England BioLabs GmbH, Hessen, Germany), and the size was determined by Bioanalyzer 2100 using Agilent High Sensitivity DNA Kit (Agilent Technologies GmbH, Waldbronn, Germany). The final library was sequenced on the Illumina MiSeq platform using v3 chemistry, generating 2 × 300 bp paired-end reads (Illumina).
Buccal DNA isolation, library preparation and sequencing of the whole mitochondrial genome
Genomic DNA isolation from buccal swab samples
Genomic DNA was isolated from the buccal swab samples using a DNeasy Blood/Tissue Kit (Qiagen) according to the standard protocol. Genomic DNA (gDNA) samples were processed for library preparation, as previously described in the Human mtDNA Genome protocol for Illumina Sequencing Platform , with small modifications. In brief, two primer sets [CLF_mtDNA_1_F (TTCTTCGGCGCATTCCACAA) and CLF_mtDNA_1_R (GGCATGCCTGTTAATGCGAG); CLF_mtDNA_2_F (TCACTTTATGCTTAGGGGCCA) and CLF_mtDNA_2_R (CAACATTTTCGGGGTATGGGC)] were used to enrich the canine mtDNA by long-range PCR. These PCR products cover the whole canine mtDNA sequencing (16,727 bp). The PCR reagents were combined in a total reaction volume of 20 µl, consisting of 10 ng template gDNA, 10 µM of each forward and reverse primer, 10 mM dNTPs, Phusion High-Fidelity DNA polymerase (Thermo Fisher, Germany) and 5× Phusion HF buffer. The PCR conditions for the amplification of the CLF_mtDNA_1 fragment were 98 °C for 30 s; 29 cycles of 98 °C for 5 s, 67.7 °C for 10 s, and 72 °C for 4 min 45 s; 72 °C for 5 min; hold at 4 °C, and those for the amplification of the CLF_mtDNA_2 fragment were 98 °C for 30 s; 29 cycles of 98 °C for 5 s, 65.8 °C for 10 s, and 72 °C for 4 min 15 s; 72 °C for 5 min; and hold at 4 °C. Library preparation was performed using a Nextera XT DNA Library Preparation Kit (Illumina Inc., CA, USA) according to the manufacturer’s instructions. The 10-pM final library was sequenced on the Illumina MiSeq sequencing platform (2 × 150 bp paired-end reads) (Illumina Inc.).
Microbiota data processing
Raw sequencing reads in fastq format were merged and filtered for low quality (vsearch  v2.12.0; maxdiffs = 2, maxee = 0.5), and chimeric sequences were removed from the data (vsearch with RDP Gold v9 as reference database). Afterwards, de-replication was performed (vsearch) and sequences were clustered into operational taxonomic units (OTUs) using UPARSE  as implemented in USEARCH (v11.0.667). Taxonomic assignment was performed applying the SINTAX  algorithm (vsearch) with RDP v18  as the classification database. Chimaera-free merged reads were mapped to OTUs using vsearch (command: usearch_global) with a 97% identity threshold to create the OTU table. The phylogenetic tree was reconstructed by aligning OTU sequences using mothur  (v.1.41.3) and SILVA  v123 as a reference database. OTUs belonging to the genus Staphylococcus were further specified at the species level using NCBI Web BLAST (Nucleotide BLAST; https://blast.ncbi.nlm.nih.gov/Blast.cgi) against the 16S rRNA database with default parameters.
Microbiome data filtering
The data were imported into R (v4.2.0) using phyloseq  (v1.40.0), and covariates were matched. OTUs not matching Bacteria or unknown phylum were excluded. Additionally, OTUs matching the phyla Abditibacteriota, Armatimonadetes, Balneolaeota, Chloroflexi, Cyanobacteria/Chloroplast, Elusimicrobia, Gemmatimonadetes, Latescibacteria, Nitrospirae, Planctomycetes, or Rhodothermaeota were excluded as potential contaminants. Next, decontam  (v1.16.0), a contamination detection tool, was applied to remove contaminants based on frequency (threshold 0.1; 114 OTUs excluded after manual inspection of potential contaminants) and prevalence (threshold 0.2; 72 OTUs excluded). Finally, only samples reaching at least 2,300 contigs were included in the downstream analysis (7 samples excluded).
Alpha and beta diversity
In faecal samples, alpha diversity was estimated using DivNet  (v0.4.0; EMiter = 6, EMburn = 3, MCiter = 250, MCburn = 100, network = diagonal). To test for significant differences in DivNet-estimated Shannon diversity, we tested for heterogeneity of total diversity (observed plus unobserved) across multiple sites using the betta  function (breakaway v4.7.9). Groupwise estimates were compared using the testDiversity function (DivNet). For skin samples, alpha diversity (Shannon index) was estimated using the estimate_richness function of phyloseq, and significance was assessed using non-parametric tests (Wilcoxon test for pairwise testing and Kruskal–Wallis test for more than 2 group comparisons with pairwise Wilcoxon tests as post hoc test).
Beta diversity was estimated using centred log-ratio transformed (clr) data and distances were calculated using Euclidean distance (Aitchison distance ). Permutational multivariate analysis of variance using distance matrices (PERMANOVA) was used to analyse differences in beta diversity (adonis2 function, as implemented in the vegan package v2.6-2, and pairwise.perm.manova function, as implemented in the RVAideMemoire package v0.9.79, each with 99,999 permutations). Inter-object dissimilarities were visualised using constrained analysis of principal coordinates (CAP) (vegan package v2.6-2).
Differentially abundant taxa
Differentially abundant (DA) taxa were identified on prevalence filtered data (taxa with a prevalence < 20% were excluded in each comparison). A count regression for correlated observations with a beta-binomial distribution model as implemented in the corncob  package (v0.2.0) was applied, and taxa with a q-value below 0.1 were considered significant. To correct for sex differences in the analysis, sex was incorporated into the NULL model in the intra-location comparisons (H0 ~ sex, H1 ~ group + sex). For the haplogroup analysis, disease group was incorporated into the NULL model to correct for these effects (H0 ~ group, H1 ~ group + haplogroup). To verify the findings, ANCOM-BC  (v1.6.0) was used (lib_cut = 0, struc_zero = TRUE, neg_lb = FALSE, tol = 1e−5, max_iter = 10000), using the same models as in the case of corncob. Again, taxa with a q value below 0.1 were considered significantly different. Differential abundance testing for Staphylococcus spp. was performed using ANCOM-BC.
If not stated differently, analyses were performed using R (v4.2.0), and p values were corrected for multiple testing using Benjamini–Hochberg correction (denoted using q or fdr). For data handling, phyloseq and tidyverse (v1.3.1) were used; cowplot (v1.1.1), ggpubr (0.4.0), sjplot (v2.8.10) and patchwork (v1.1.1) were applied to create the figures.
Balances with CADESI-4 or PVAS as the target variable were calculated using a forward-selection method as implemented in selbal  (v0.1.0) on phyla/genera occurring in at least one-third of the samples per location with 5-fold cross-validation (10 iterations each).
For correlation analysis (Spearman rank correlation), clr-transformed abundance values were used. At the genus level, only genera found in at least 1/15th of the samples were included to calculate correlations.
The workflow used for the microbiota analysis is available at https://github.com/buschlab/2022_canine_atopic_dermatitis_paper.
Mitochondrial genome analysis
Raw sequencing data were mapped with bwa (version 0.7.17)  against dog reference sequence U96639.2 (https://www.ncbi.nlm.nih.gov/nuccore/U96639), and variant calling was performed with bcftools (version 1.15)  using parameter --ploidy 1, i.e., assuming a ploidy of 1. The PHRED-scaled genotype likelihood was visualised using a heatmap. This value ranges between 0 and 255, with 0 being the strongest support of a homomorphic non-reference (alternative) allele. Values larger than 0 were typically seen at positions at which more than one allele is observed, indicative of variant calling artefacts or heteroplasmic variants. Variants with no genotype calls, i.e., homomorphic reference allele, are displayed in white. Mitochondrial (MT) sequences were generated by removing indels and variants with quality less than 200 and then applying the remaining variants to the reference sequence using the bcftools consensus. The resulting MT sequences have individually been uploaded to the Canis mtDNA HV1 database (http://chd.vnbiology.com)  for haplogroup assignment. The MT analyses workflow used is available at https://github.com/TLenfers/multispecies_mitochondrial_variant_analysis. Furthermore, MT sequences used for haplogroup assignment as well as the code for generating Figure S11, including underlying variant data, are available at https://github.com/iwohlers/2022_dog_mt.
Differential microbiota composition in healthy and cAD-affected dogs with and without medical treatment
Five hundred fifty-four skin swab samples were successfully sequenced for the hypervariable V1-V2 region of the 16S rRNA gene. On average, 44,387 (s.d. ± 26,307) read pairs were obtained in each sample. Merging forward and reverse reads into contigs, quality filtering, and chimaera removal resulted in an average of 29,343 (± 19,755) contigs per sample.
Taxonomic classification of 12 skin sites from healthy dogs, untreated cAD-affected dogs (cADpre), the same cADpre dogs that received 2 weeks of oclacitinib treatment (cADpost), and cAD-affected dogs that had already been on the treatment (cADtreat) dogs at the phylum and genus levels are shown in Fig. 1a and Fig. 1b (bacterial taxa identified in these skin samples are shown in Supplementary Material S2 as well). While there is an obvious heterogeneity in the relative abundance of the skin microbiota across the 12 skin sites, Staphylococcus was the predominant bacterial taxon in all 12 skin sites of cAD-affected dogs (Fig. 1b). More specifically, its abundance in cADpre dogs was significantly higher than that in healthy dogs (Healthy) in abdomen (p = 0.0169), axilla (p = 0.0330), cubital flexor (p < 0.0001), palmar metacarpal (p = 0.0157), and perilabial area (p = 0.0153). In abdomen, cubital flexor, perilabial area, 2 weeks of treatment with JAK inhibitor significantly reduced the abundance of Staphylococcus (p = 0.0419, p = 0.0005, p = 0.0299, respectively; cADpre vs cADpost; Fig. 1c).
Reduction of bacterial diversity in cAD skin
Skin microbiome data of the 12 skin sites were analysed for alpha diversity (Figure S2), and the difference in skin bacterial community composition was evaluated by beta diversity (Figure S3). Alpha diversity (estimated using Shannon diversity) showed a general trend of reduction of bacterial richness in untreated cAD-affected dogs (cADpre) compared with healthy dogs, and the level was recovered after the medical intervention. Of the 12 skin sites, a significant difference between the four groups was observed in axilla (p = 0.01), medial pinnae (p = 0.023) and perilabial area (p = 0.023). A significant reduction of alpha diversity in cAD-affected dogs (cADpre) compared with healthy dogs (Healthy) was observed in perilabial area (p = 0.011), cubital flexor (p = 0.041) and medial pinnae (p = 0.006). When the same dogs before (cADpre) and after the 2-week JAK inhibitor treatment (cADpost) was compared, a significant restoration of alpha diversity was observed in medial pinnae (p = 0.019). After the longer term treatment of JAK inhibitor, the alpha diversity
Beta diversity was significantly different between groups in abdomen, front paw, perineum, and ventral tail (Figure S3).
Differentially abundant bacterial taxa in the skin between healthy and cAD-affected dogs with and without oclacitinib treatment
For further analysis, we selected the skin sites that were most relevant to the clinical score. To select these skin sites, balances of bacterial taxa were calculated using a greedy stepwise algorithm (selbal) with 5-fold cross-validation and 10 iterations. This method estimates the optimal number of discriminating taxonomic groups and the ratio of these taxonomic groups (named “denominator” and “numerator”) and then defines the balance between the microbial characteristics that best describe the differences in the compared phenotype—in this case, CADESI-04. The skin sites whose bacterial taxa at the genus level significantly correlated with their CADESI-04 scores (R2 > 0.4) were selected for further analysis (Table S2). To this end, skin microbiome data from four skin sites, namely, the abdomen, front paw, perilabial area, and ventral tail, were selected for more detailed analysis.
For the four selected skin sites, differentially abundant bacterial taxa at the phylum level were analysed between two groups: cADpre vs. healthy, cADtreatment vs. healthy, and cADpost vs. cADpre (Figure S4 and Table S3). The abundance of Firmicutes in the front paw and perilabial area was significantly higher in cAD-affected dogs without treatment than in healthy dogs (cADpre vs. healthy). At these skin sites, cAD-affected dogs showed a significant reduction in Firmicutes abundance after 2 weeks of oclasitinib treatment (cADpost vs. cADpre). At all four skin sites, the abundance of Proteobacteria was significantly increased in cAD-affected dogs after they received 2 weeks of oclacitinib treatment (cADpost vs. cADpre). In contrast, the abundance of Deinococcus-Thermus significantly decreased after treatment in the abdomen, front paw, and ventral tail. Similarly, differentially abundant bacterial taxa at the genus level were analysed between the same two groups (Figure S5 and Table S4). The abundance of Staphylococcus was significantly higher in the perilabial area and ventral tail of cAD-affected dogs than in the same skin sites of healthy dogs. When compared with paired samples before and after the 2-week oclacitinib treatment (cADpost vs. cADpre), the abundance of Staphylococcus in the perilabial area was largely decreased after the treatment (q < 1.00 × 10−50, effect size = − 2.4017), while the change in abdomen, front paw, and ventral tail was marginal (q = 3.50 × 10−41, effect size = 0.8125; q = 1.64 × 10−6, effect size = 0.1715; q = 2.25 × 10−8, effect size = − 0.3779, respectively). Interestingly, the abundance of Staphylococcus on the abdomen of cAD-affected dogs that already received oclacitinib treatment (cADtreat), was less than that in the site-matched skin of healthy dogs. Meanwhile, there was still a higher bacterial load on the ventral tail of cADtreat dogs than in the same skin site on healthy dogs (q = 0.0554, effect size = 1.5864).
Correlation analysis between bacterial taxa and clinical cAD parameters
Next, we conducted a correlation analysis between clinical parameters (CADESI-04 and PVAS) and microorganisms on the skin in each group to evaluate whether identified differences in the abundance of bacteria were associated with the clinical parameters (Fig. 2 and Figure S6, respectively). To identify specific bacterial phyla and/or genera that correlate with respective clinical phenotypes and are commonly shared in all cAD-affected dogs, we looked for similar trends of correlation in all three dog groups in all four skin sites. More specifically, bacterial taxa showing similar colours across the three groups in each heatmap were evaluated. There were no bacterial phyla and/or genera that had a similar impact (i.e., negative or positive relations) in CADESI-04 or PVAS. Only two genera, Corynebacterium and Staphylococcus, demonstrated a trend of positive correlation with CADESI-04 values in all groups in all four skin sites (Fig. 2b).
Association of Staphylococcus abundance with cAD
Having our correlation analysis as well as previously reported clinical evidence, we further investigated whether the abundance of Staphylococcus was associated with cAD in all 12 skin sites (Figure S7 and Table S5). Of these skin sites, five (namely, the abdomen, axilla, cubital flexor, palmar metacarpal, and perilabial area) showed a significant association of Staphylococcus abundance between cADpre dogs and healthy dogs (p = 0.0392, p = 0.0436, p = 0.0026, p = 0.0105, and p = 0.0233, respectively; Figure S7a–e). Even though the p values were not statistically significant (due to the large variation between individuals), the abundance of Staphylococcus in the other three skin sites (hind paw, flank, and front paw) showed a trend towards higher abundances in cADpre dogs than in healthy dogs (p = 0.2331, p = 0.1221, and p = 0.2889, respectively; Figure S7f–h). These skin sites demonstrated a positive correlation with the corresponding CADESI-04 value (Figure S7i–p).
Of the identified OTUs, 14 different OTUs were mapped to the genus Staphylococcus. These OTUs were further classified at species level by an NCBI Nucleotide BLAST search. The ten most abundant species in 12 skin sites are presented in Figure S8 and Table S6. In 9 of the 12 skin sites, the abundance of S. pseudintermedius was significantly higher in cAD-affected dogs not receiving medical treatment (cADpre) than in healthy dogs (q < 0.1; Fig. 3 and Table S7).
Distinct gut microbial composition between healthy and cAD-affected dogs with and without oclacitinib treatment
A total of 46 faecal samples were analysed. Approximately 41,000 (± 25,895) read pairs per sample were obtained, and merging forward and reverse reads into contigs yielded approximately 41,000 (± 21,305) contigs per sample. After filtering and chimaera removal, the average number of contigs per sample was 36,367 (± 17,293).
The bacterial taxonomic composition in the faecal samples from all dogs was plotted (Fig. 4). The mean relative abundances of bacterial taxa in the four groups were investigated. A total of ten phyla were identified in the faecal samples (Fig. 4a). Firmicutes was the most abundant, followed by Bacteroidetes and Proteobacteria. We observed a major decrease in the proportion of Fusobacteria accompanied by an increase in Proteobacteria in all cAD groups (cADpre, cADpost, and cADtreat) (Table S8).
A total of 131 different genera were identified in the faecal samples. The mean relative abundance of the 19 most abundant genera is shown in Fig. 4b. In healthy dogs, the predominant genera were Fusobacterium (20.6%), Megamonas (18.4%), and an unclassified genus belonging to the Bacteroidaceae family (16.2%). In contrast, the predominant genera found in untreated cAD-affected dogs (cADpre) were Escherichia/Shigella (19.4%), Bacteroides (18.1%), and Clostridium sensu stricto (9.7%), with Fusobacterium (0.06%) and Megamonas (0.0003%) being detected at very low abundances in this group (Table S9).
Next, differentially abundant bacterial taxa in the faecal samples at the phylum level were analysed between pairs of groups; healthy vs. cADpre, cADtreat vs. healthy, and cADpre vs. cADpost (Figure S9, Table 1). In untreated cAD-affected (cADpre) dogs, the abundance of Fusobacteria was significantly lower than that in healthy dogs (effect size = − 5.3705, q = 1.38 × 10−16), while the abundance was increased after 2 weeks of oclacitinib treatment (effect size = 5.11224, q < 0.0001). No change in this specific bacterial phylum was observed when cAD-affected dogs that already received oclacitinib treatment (cADtreat) and healthy dogs were compared.
Similarly, differentially abundant faecal bacterial taxa between the groups at the genus level were analysed (Fig. 4c and Table 2). The abundance of Fusobacterium was significantly lower in untreated cAD (cADpre) dogs than in healthy dogs (effect size = − 5.3391, q = 3.64 × 10−10). It was increased in the same cAD-affected dogs after they received 2-week-oclacitinib treatment (cADpost vs. cADpre; effect size = 5.1177, q < 0.0001). A similar trend (increase in cADpre and decrease after the 2-week treatment) was observed in Blautia (effect size = −1.1738, q = 3.47 × 10−245), Clostridium sensu stricto (effect size = − 1.2619, q = 5.67 × 10−58) and Escherichia/Shigella (effect size = −0.4567, q = 8.7837 × 10−106). Dogs in the cADtreat group did not show a reduction in the abundance of Blautia and Clostridium sensu stricto, which remained higher than that in healthy dogs (Blautia, q = n.s; C. sensu stricto, effect size = 1.6453, q = 4.80 × 10−28).
Correlation between faecal bacterial taxa and clinical cAD parameters
Next, we conducted a correlation analysis between clinical parameters (CADESI-04 total score and PVAS; Figure S10) and bacterial taxa in the gut in each group. To evaluate specific bacterial phyla and/or genera correlating with respective clinical phenotypes commonly shared in all dogs, we looked for similar trends of correlation in all four dog groups. No bacterial phyla or other taxa were observed to correlate with PVAS (Figure S10a, b). Three phyla (Bacteroidetes, Actinobacteria, and Firmicutes) and three genera (Amedibacillus, Bacteroides, and Flavonifractor) were positively correlated with CADESI-04 total score, while three phyla (Acidobacteria, Proteobacteria, and Synergistetes) and four genera (Cellulosilyticum, Escherichia/Shigella, Novosphingobium, and Stenotrophobacter) were inversely correlated with the score (Figure S10c, d).
Differential gut microbial diversity in cAD-affected dogs compared with healthy dogs was reversed by oclacitinib treatment
Species diversity in the faecal samples was calculated using the DivNet estimate of Shannon that integrates the number of different species present (species richness) as well as their relative distribution (species evenness). DivNet also accounts for interactions between bacterial communities in the ecological network and for missing taxa. When alpha diversity was measured samplewise, no significant differences were observed between any of the four groups (Fig. 4d). Samples from all groups showed a wide range of species diversities inside their group, with mean values of 2.01 (healthy), 1.88 (cADpre), 1.96 (cADpost), and 2.06 (cADtreat). Therefore, there was a tendency towards a reduction in species diversity in untreated cAD-affected dogs (cADpre) compared with healthy dogs. Next, we conducted a communitywise estimate, which resulted in statistically significant differences between the groups (Fig. 4e). In this analysis, all samples per group are combined and regarded as one community so that the overall species diversity of the healthy or diseased environment is estimated. Compared with healthy dogs, community diversity was significantly decreased in untreated cAD (cADpre) dogs (healthy, 2.06 ± 0.35, p < 0.001; ADpre, 1.01 ± 0.73, p < 0.001). Similarly, the community after 2 weeks of treatment was less diverse compared to the healthy group (0.96 ± 0.57, p < 0.001). In the cAD groups that already received oclacitinib treatment (cADtreat), the mean diversity became closer to that of healthy dogs (1.92 ± 0.66), although a larger deviation between individuals was observed.
Beta diversity measures dissimilarity between the bacterial composition of multiple populations. The “Aitchison distance”, which is the Euclidian distance after clr (centred log-ratio) transformation, was used to measure beta diversity and constrained analysis of principal coordinates to visualise and evaluate clustering of the samples according to study group and disease status (Fig. 4f).
To analyse the ability of the samples to be separated by the group (Healthy, cADpre, cADpost and cADtreat) and disease severity (CADESI-04 score), multivariate analysis of variance using the distance matrices (adonis) was used. The disease separate categories were applied according to the published consensus , i.e., none, < 10; mild cAD, 10–34; moderate cAD, 35–59, and severe cAD, ≥ 60. Aitchison distances showed that the group explained 14.6% (p < 0.0001) and disease severity explained 11.5% (p < 0.001) of the variance (r2) and that the healthy dogs were found to be significantly different compared with each of the other cAD groups (PERMANOVA, p < 0.001 for cADpre and cADpost; p < 0.001 for cADtreat). When testing for disease severity, beta diversity significantly differed between disease-free animals (none) and those with mild cAD (p < 0.05), while other comparisons did not reach statistical significance. Neither the sex (adonis, p = 0.64) of the dogs nor the age group (adonis, p = 0.09) was found to be significant in beta diversity.
In summary, the bacterial community in all dogs with cAD, regardless of the treatment status, was significantly different from that in healthy dogs.
Mitochondrial DNA variant clustering and mitochondrial DNA haplogroup determination
Buccal swab samples were obtained, and gDNA was isolated from all 40 dogs (20 healthy, 10 cADpre, and 10 cADtreat). Two swab samples did not yield sufficient gDNA for further processing. Ultimately, samples from 20 healthy dogs, nine cADpre dogs, and nine cADtreatment dogs were successfully sequenced for whole canine mitochondrial DNA. The obtained sequencing data were aligned to the reference sequence and variant calling was performed.
Figure S11 shows the complete mitochondrial genetic profiles of all 38 dogs. For haplogroup determination, we first used the Canis mtDNA HV1 database, which is available online (http://chd.vnbiology.com). The haplogroups determined by this approach overlapped the heatmap clustering pattern. Therefore, we decided to use major haplogroup A (combined all A sub-haplogroups, except A64; a total of 20 dogs, of which 16 had microbiome data available) and haplogroup C17 (including three samples assigned “new haplogroup C ref: C17”, thus together referred to as C afterwards; a total of 14 dogs, of which 12 dogs had microbiota data available) for further association analysis.
Mitochondrial DNA haplogroups are associated with the gut and skin microbiome
We hypothesised that genetic variation in the mtDNA is associated with cAD. Since all homoplasmic variants in the tested dogs are haplogroup-specific (Figure S11), we performed a χ2 test for the mtDNA haplogroup distribution between healthy and diseased samples (cADpre and cADtreat). There was no association between the mtDNA haplogroup and cAD in this cohort (p = 0.2675).
We further investigated the potential association between mitochondrial haplogroups and microbiota in the gut and the skin.
We first conducted the association between mtDNA haplogroup and gut microbiota. To evaluate the association between mtDNA haplogroup and alpha diversity, a random effects model (no random slope/intercept) was used with group (healthy, cADpost, and cADtreat) as a random effect. The model shows that the effect of the haplogroup on alpha diversity is not significant (p = 0.1194; Fig. 5a). Additionally, testing for heterogeneity of total diversity revealed a significant reduction in alpha diversity in haplogroup C (Fig. 5b; beta: p = 0.033, estimate = − 0.3730).
Beta diversity was estimated using Aitchison distance, and a constrained correspondence analysis was applied to estimate differences and correct for group effects (estimating the haplogroup effects; Fig. 5c). The model was not significant (p = 0.07822). Neither haplogroup (p = 0.07845) nor the separation at the x-axis (CAP1, p = 0.0784) were significant; p values were calculated using 99,999 permutations using an ANOVA-like permutation test for constrained correspondence analysis. Although these values were not statistically significant, they did show a trend (p < 0.1).
Finally, we analysed the differential abundance using corncob analysis (Fig. 5d, e). Three phyla (Actinobacteria, Campylobacterota, and Firmicutes) and seven genera (Amedibacterium, Anaerostipes, Anaerotignum, Blautia, Clostridium XVII, Roseburia, and Streptococcus) were significantly lower in the C haplogroup than in the A haplogroup (q < 0.1; Table S10). For this analysis, group was included in the null model to remove group effects from the analysis.
The same approach was taken to evaluate the skin microbiota and mtDNA haplogroups. For the skin samples, we selected and evaluated three skin sites, namely axilla (Figure S12), medial pinnae (Figure S13) and front paw (Figure S14), which exhibited a significant difference in alpha diversity between the four dog groups. No significant correlation between mtDNA haplogroups and alpha diversity (Figure S12a,b, Figure S13a,b and Figure S14a,b) or beta diversity (Figure S12c, Figure S13c and Figure S14c) was observed in any of the three skin sites. When the differential abundance of bacteria was analysed, we observed the significant difference between the haplogroups (Figures S12d,e, S13d,e and S14d,e; Table S10); In axilla, two phyla (Actinobacteria and Bacteroides) and four genera (Arthrobacter, Bacteroides, Rubellimicrobium and Haemophilus) were significantly differentially abundant between C and A haplogroups (q < 0.1). In medial pinnae, one phylum Actinobacteria was significantly reduced in haplogroup A than haplogroup C (q < 0.1). At the genus level, Flavobaterium and Pasteurella were higher in haplogroup C than haplogroup A (q < 0.1), while Rothia, Stenotrophomonas and Tannerella were higher in haplogroup A than haplogroup C (q < 0.1). In front paw, four genera (Bacteroides, Capnocytophaga, Neisseria, and Haemophilus) were differentially abundant between the two mtDNA haplogroups (q < 0.01).
In summary, there was a trend of differences in beta diversity in faecal samples, and we observed significantly different abundant phyla/genera both in faeces and skin samples between mtDNA haplogroups.
The skin and commensal microbiota are critical players in the health and pathological status of hosts. AD is an allergic skin disorder observed in both humans and dogs.
In this study, 16S rRNA gene amplicon sequencing of skin swabs and faecal samples obtained from dogs affected with cAD and healthy Shiba Inu dogs was conducted.
From the skin microbiota analysis, we demonstrated clear dysbiosis on the skin of cAD-affected dogs, with a higher abundance of Staphylococcus in all 12 skin sites of cAD-affected Shiba Inu dogs. These findings are in line with previously reported studies. The paired samples, i.e., untreated cAD (cADpre) and the same cAD-affected dogs after 2 weeks of oclacitinib treatment (cADpost), demonstrated that Staphylococcus abundance, predominantly S. pseudintermedius, was significantly lower in cADpost dog skin than in cADpre dog skin. Upon oclacitinib treatment, the abundance was significantly reduced. At the same time, we conducted the analysis of faecal microbiota in healthy dogs and cAD-affected dogs, in parallel to the analysis of their skin microbiota. Additionally, we were able to evaluate the therapeutic impact of oclacitinib on the gut microbiota of untreated cAD-affected dogs that received 2 weeks or cAD dogs that were already on the treatment for an average of 13 months with oclacitinib. In fact, because untreated cAD-affected dogs (cADpre) and 2-week oclacitinib-treated cAD-affected dogs were paired, we were able to evaluate the gut microbiota in pre- and post- treatment cAD status.
In parallel to the skin microbiota, there was clear dysbiosis in the gut of cAD-affected dogs. In healthy dogs, the predominant phyla in the faecal samples included Firmicutes, Bacteroidetes, Proteobacteria and Fusobacteria. This observation is in line with previous studies . When comparing the gut microbiota in healthy dogs with that in cAD-affected dogs (cADpre), the abundance of Fusobacteria was significantly reduced in cAD-affected dogs, which is in line with previous studies of gut microbiota in dogs affected with cAD and healthy dogs [22, 23]. This was further illustrated at the genus level, presenting the reduction of Fusobacterium abundance in cADpre dogs compared with that in healthy dogs. In fact, Fusobacterium was primarily detected in several carnivores (African lion, arctic fox, king vulture) and non-carnivores (dogs and meerkat) as commensal microbiota . In humans, Fusobacterium is known to be highly abundant in various gastrointestinal diseases, such as inflammatory bowel disease (IBD) , ulcerative colitis , and colorectal cancer , suggesting that the impact of Fusobacterium on health and disease in dogs seems to be different from that in humans . Megamonas, which belongs to the phylum Firmicutes (Bachillota), is another bacterial genus that was significantly less abundant in cADpre dogs than in healthy dogs. This finding is also consistent with the previous cAD study conducted in the same area as our study . Another pilot study of gut microbiota in cAD dogs showed the opposite finding; a higher prevalence of Megamonas in faecal samples from four cAD-affected beagle dogs compared with three healthy beagle dogs . As such, differential findings between studies could be due to the environmental factors (e.g., geographical location and diet), dog breeds as well as sample sizes used in the studies. These bacterial taxa are known to be pathologically relevant, such as in IBD in dogs; more specifically, the abundances of both Fusobacterium and Megamonas were negatively correlated with canine IBD . In dogs, the abundance of Fusobacterium is positively correlated with faecal concentrations of butyrate and propionate , and some species of Megamonas can ferment glucose into acetate and propionate [60, 61]. Given that the higher levels of butyrate and propionate in faeces reduced the risk of atopy in early life in humans, the reduction of such SCFA-producing bacteria in cAD-affected dogs suggests that Fusobacteria and Megamonas could explain the pathological relevance of cAD in Shiba Inu dogs. Furthermore, several meta-analysis studies of the published data suggested an association between IBD and AD in humans, i.e., a higher incidence of AD in IBD patients as well as a higher prevalence of IBD in AD patients [62, 63], suggesting that the phylum Fusobacteria and the genus Fusobacterium may be potential biomarkers not only for the risk of cAD but also for the risk of IBD in dogs. Nevertheless, further studies of faecal microbiota from many different dog breeds with larger sample sizes will be needed to confirm this hypothesis.
In contrast, the abundances of Clostridium sensu stricto (former Clostridium cluster I) and Escherichia/Shigella were significantly increased in untreated cAD-affected dogs (cADpre) compared with healthy dogs. Clostridium was reported that its colonisation in the gut at an early age (at ages 5 and 13 weeks) was associated with a higher risk of AD in children  and that the higher abundance of Clostridium spp. is associated with AD in infants. Two bacterial taxa, Escherichia coli and Shigella spp., were not differentiated by the 16S rRNA gene sequence because of their genetic relatedness [65, 66]. Despite some differential biochemical characteristics, both bacterial taxa are enteric pathogens, and their virulence factors influence many cellular processes . Higher abundances of both Shigella and E. coli in faeces were reported to be associated with atopic eczema [67, 68].
These changes observed in the gut microbiota of the untreated cAD-affected dogs (cADpre) were further altered most likely by medical intervention with oclacitinib. After the 2-week treatment, the abundance of Fusobacterium, which was less than 2% in cAD-affected dogs, was significantly increased in the same dogs (cADpost vs. cADpre). On the other hand, the abundance of Escherichia/Shigella was significantly decreased after the 2-week treatment. In cAD-affected dogs that were already on oclacitinib treatment (cADtreat), the composition of the gut bacterial taxa was richer than that of the cAD-affected dogs that received 2-week treatment (cADpost): Fusobacterium was retained, the abundance of Escherichia/Shigella was reduced, and Megamonas appeared. This finding indicates that the 2-week oclacitinib treatment shifted gut microbial composition, including the increased abundance of Fusobacterium, which may be beneficial for dog health. The gut microbial composition of the faeces from cAD-affected dogs that were already on oclacitinib treatment (cADtreat) showed relatively similar patterns compared with untreated cAD-affected dogs (cADpre), i.e., higher abundances of Escherichia/Shigella and Clostridium sensu stricto, yet we observed higher abundances of Fusobacterium and Megamonas, which were also highly abundant in healthy dogs, suggesting that the treatment with oclacitinib may shift the gut microbial composition of cAD-affected dogs towards that of healthy dogs to a certain extent. Five out of 10 cADtreat dogs that already received oclacitinib treatment, were on a varied prescription diet specifically designed for allergy control; thus, this dietary control may have also supported modulating the gut microbiota composition from “cAD gut type” to “healthy gut type”.
We also sequenced the whole mitochondrial genome of Shiba Inu dogs with and without cAD in this study. In humans, mitochondrial haplogroups are utilised in the context of population genetics because they represent a phylogenetic tree that allows retracing human settlement of the world, resulting in present-day worldwide geographic differences in haplogroup abundances. In contrast, dogs have been recently crossed, and their breeding has been largely controlled by humans. Therefore, we hypothesised that mtDNA variants may be variable among one specific breed, which may contribute to host gene-microbiota interactions in the pathology of cAD. We did not observe a significant association between mtDNA haplogroups and cAD in this sample cohort; however, this may be due to the limited sample size. Yet, we identified a significant association between mtDNA haplogroups and specific bacterial taxa. Further studies using larger sample sizes may indicate that the mtDNA haplogroup and specific gut microbes are biomarkers to evaluate disease predisposition, treatment efficacy, and health management in dogs.
In our present study, bacterial 16S rRNA gene sequencing was used for both skin and faecal microbiome analysis. Despite its low cost and availability of well-established analysis pipelines, this approach has limitations in taxonomic resolution and low sensitivity in comparison to whole genome shotgun metagenome sequencing . For example, this study did not cover the fungal community in the skin and faeces in dogs. Colonisation of fungi, especially, Malassezia species in the skin, is well-known to be associated with cAD . Detailed identification of bacterial species and functional profiling of identified microbes are of interest, especially Staphylococcus species. Therefore, in our extended future study, we will conduct shotgun metagenomics sequencing to cover these aspects in both skin and faecal samples from healthy and cAD-affected dogs.
We conducted a comprehensive microbial analysis in the skin and gut of cAD and healthy Shiba Inu dogs. The cAD-affected dogs included untreated cAD-affected dogs, the same dogs after a 2-week course of oclacitinib treatment, and an independent group of cAD-affected dogs that had already been receiving oclacitinib treatment for up to two years. This study is the first report to present clear dysbiosis of both the skin and the gut in dogs with cAD. We observed that oclacitinib treatment on cAD, regardless of the treatment duration, shifted both the skin and gut microbiota toward those of healthy dogs. In addition, we identified a significant association between mtDNA haplogroups and specific bacterial taxa in the gut and the skin microbiota in dogs. This finding will be the basis for novel disease management strategies for cAD.
Availability of data and materials
Skin and faecal 16S rRNA sequencing data and dog mitochondrial genome sequencing data were submitted to the European Nucleotide Archive (ENA) and are available under accession numbers PRJEB53048 (microbiome data) and PRJEB53275 (whole mitochondrial genome sequencing data).
Canine atopic dermatitis
Canine atopic dermatitis evaluation score index-04
Pruritis visual analogue scale
16S ribosomal RNA
Hayward JJ, Castelhano MG, Oliveira KC, Corey E, Balkman C, Baxter TL, et al. Complex disease and phenotype mapping in the domestic dog. Nat Commun. 2016;7:10460.
LeBlanc AK, Mazcko CN. Improving human cancer therapy through the evaluation of pet dogs. Nat Rev Cancer. 2020;20:727–42.
Ruple A, MacLean E, Snyder-Mackler N, Creevy KE, Promislow D. Dog models of aging. Annu Rev Anim Biosci. 2022;10:419–39.
Weidinger S, Beck LA, Bieber T, Kabashima K, Irvine AD. Atopic dermatitis. Nat Rev Dis Primer. 2018;4:1–20.
Gedon NKY, Mueller RS. Atopic dermatitis in cats and dogs: a difficult disease for animals and owners. Clin Transl Allergy. 2018;8:41.
Nuttall TJ, Marsella R, Rosenbaum MR, Gonzales AJ, Fadok VA. Update on pathogenesis, diagnosis, and treatment of atopic dermatitis in dogs. J Am Vet Med Assoc. 2019;254:1291–300.
Marsella R, Girolomoni G. Canine models of atopic dermatitis: a useful tool with untapped potential. J Invest Dermatol. 2009;129:2351–7.
Jaeger K, Linek M, Power HT, Bettenay SV, Zabel S, Rosychuk RAW, et al. Breed and site predispositions of dogs with atopic dermatitis: a comparison of five locations in three continents. Vet Dermatol. 2010;21:118–22.
Williams MR, Gallo RL. The role of the skin microbiome in atopic dermatitis. Curr Allergy Asthma Rep. 2015;15:65.
Carding S, Verbeke K, Vipond DT, Corfe BM, Owen LJ. Dysbiosis of the gut microbiota in disease. Microb Ecol Health Dis. 2015;26:26191.
Kong HH, Oh J, Deming C, Conlan S, Grice EA, Beatson MA, et al. Temporal shifts in the skin microbiome associated with disease flares and treatment in children with atopic dermatitis. Genome Res. 2012;22:850–9.
Rodrigues Hoffmann A, Patterson AP, Diesel A, Lawhon SD, Ly HJ, Elkins Stephenson C, et al. The skin microbiome in healthy and allergic dogs. PloS One. 2014;9:e83197.
Bradley CW, Morris DO, Rankin SC, Cain CL, Misic AM, Houser T, et al. Longitudinal evaluation of the skin microbiome and association with microenvironment and treatment in canine atopic dermatitis. J Invest Dermatol. 2016;136:1182–90.
Chermprapai S, Ederveen THA, Broere F, Broens EM, Schlotter YM, van Schalkwijk S, et al. The bacterial and fungal microbiome of the skin of healthy dogs and dogs with atopic dermatitis and the impact of topical antimicrobial therapy, an exploratory study. Vet Microbiol. 2019;229:90–9.
Fazakerley J, Nuttall T, Sales D, Schmidt V, Carter SD, Hart CA, et al. Staphylococcal colonization of mucosal and lesional skin sites in atopic and healthy dogs. Vet Dermatol. 2009;20:179–84.
Lee SY, Lee E, Park YM, Hong SJ. Microbiome in the gut-skin axis in atopic dermatitis. Allergy Asthma Immunol Res. 2018;10:354–62.
O’Neill CA, Monteleone G, McLaughlin JT, Paus R. The gut-skin axis in health and disease: A paradigm with therapeutic implications. BioEssays News Rev Mol Cell Dev Biol. 2016;38:1167–76.
Smith PM, Howitt MR, Panikov N, Michaud M, Gallini CA, Bohlooly-Y M, et al. The microbial metabolites, short-chain fatty acids, regulate colonic Treg cell homeostasis. Science. 2013;341:569–73.
Schwarz A, Bruhs A, Schwarz T. The short-chain fatty acid sodium butyrate functions as a regulator of the skin immune system. J Invest Dermatol. 2017;137:855–64.
Fang Z, Li L, Zhang H, Zhao J, Lu W, Chen W. Gut microbiota, probiotics, and their interactions in prevention and treatment of atopic dermatitis: a review. Front Immunol. 2021;12:720393.
Rostaher A, Morsy Y, Favrot C, Unterer S, Schnyder M, Scharl M, et al. Comparison of the gut microbiome between atopic and healthy dogs-preliminary data. Anim Open Access J MDPI. 2022;12:2377.
Uchiyama J, Osumi T, Mizukami K, Fukuyama T, Shima A, Unno A, et al. Characterization of the oral and faecal microbiota associated with atopic dermatitis in dogs selected from a purebred Shiba Inu colony. Lett Appl Microbiol. 2022;75:1607–16.
Sugita K, Shima A, Takahashi K, Ishihara G, Kawano K, Ohmori K. Pilot evaluation of a single oral fecal microbiota transplantation for canine atopic dermatitis. Sci Rep. 2023;13:8824.
Tanaka K, Yamamoto-Fukuda M, Takizawa T, Shimakura H, Sakaguchi M. Association analysis of non-synonymous polymorphisms of interleukin-4 receptor-α and interleukin-13 genes in canine atopic dermatitis. J Vet Med Sci. 2020;82:1253–9.
Fujimura M, Ishimaru H, Nakatsuji Y. Fluoxetine (SSRI) treatment of canine atopic dermatitis: a randomized, double-blind, placebo-controlled, crossover trial. Pol J Vet Sci. 2014;17:371–3.
Wood SH, Ke X, Nuttall T, McEwan N, Ollier WE, Carter SD. Genome-wide association analysis of canine atopic dermatitis and identification of disease related SNPs. Immunogenetics. 2009;61:765–72.
Ma J, Coarfa C, Qin X, Bonnen PE, Milosavljevic A, Versalovic J, et al. mtDNA haplogroup and single nucleotide polymorphisms structure human microbiome communities. BMC Genomics. 2014;15:257.
Hirose M, Künstner A, Schilf P, Sünderhauf A, Rupp J, Jöhren O, et al. Mitochondrial gene polymorphism is associated with gut microbial communities in mice. Sci Rep. 2017;7:15293.
Hensel P, Santoro D, Favrot C, Hill P, Griffin C. Canine atopic dermatitis: detailed guidelines for diagnosis and allergen identification. BMC Vet Res. 2015;11:196.
Olivry T, Marsella R, Iwasaki T, Mueller R, Dermatitis TITFOCA. Validation of CADESI-03, a severity scale for clinical trials enrolling dogs with atopic dermatitis. Vet Dermatol. 2007;18:78–86.
Olivry T, Saridomichelakis M, Nuttall T, Bensignor E, Griffin CE, Hill PB, et al. Validation of the Canine atopic Dermatitis Extent and Severity Index (CADESI)-4, a simplified severity scale for assessing skin lesions of atopic dermatitis in dogs. Vet Dermatol. 2014;25(77–85):e25.
Reimer-Taschenbrecker A, Künstner A, Hirose M, Hübner S, Gewert S, Ibrahim S, et al. Predominance of staphylococcus correlates with wound burden and disease activity in dystrophic epidermolysis bullosa: a prospective case-control study. J Invest Dermatol. 2022;142:2117-2127.e8.
Künstner A, Schilf P, Busch H, Ibrahim SM, Hirose M. Changes of gut microbiota by natural mtDNA variant differences augment susceptibility to metabolic disease and ageing. Int J Mol Sci. 2022;23:1056.
Human mtDNA Genome For the Illumina Sequencing Platform. Available online at https://emea.support.illumina.com/content/dam/illumina-support/documents/documentation/chemistry_documentation/samplepreps_legacy/human-mtdna-genome-guide-15037958-01.pdf. Accessed 21 Sep 2023.
Rognes T, Flouri T, Nichols B, Quince C, Mahé F. VSEARCH: a versatile open source tool for metagenomics. PeerJ. 2016;4:e2584.
Edgar RC. UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat Methods. 2013;10:996–8.
Edgar RC. SINTAX: a simple non-Bayesian taxonomy classifier for 16S and ITS sequences. bioRxiv. 2016. https://doi.org/10.1101/074161.
Cole JR, Wang Q, Fish JA, Chai B, McGarrell DM, Sun Y, et al. Ribosomal Database Project: data and tools for high throughput rRNA analysis. Nucleic Acids Res. 2014;42 Database issue:D633-642.
Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol. 2009;75:7537–41.
Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41:D590-6.
McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PloS One. 2013;8:e61217.
Davis NM, Proctor DM, Holmes SP, Relman DA, Callahan BJ. Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. Microbiome. 2018;6:226.
Willis AD, Martin BD. Estimating diversity in networked ecological communities. Biostatistics. 2020. https://doi.org/10.1093/biostatistics/kxaa015.
Willis A, Bunge J, Whitman T. Improved detection of changes in species richness in high diversity microbial communities. J R Stat Soc Ser C Appl Stat. 2017;66:963–77.
Aitchison J. The statistical analysis of compositional data. J R Stat Soc Ser B Methodol. 1982;44:139–77.
Martin BD, Witten D, Willis AD. Modeling microbial abundances and dysbiosis with beta-binomial regression. Ann Appl Stat. 2020;14:94–115.
Lin H, Peddada SD. Analysis of compositions of microbiomes with bias correction. Nat Commun. 2020;11:3514.
Susin A, Wang Y, Lê Cao K-A, Calle ML. Variable selection in microbiome compositional data analysis. NAR Genomics Bioinforma. 2020;2:lqaa029.
Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinforma Oxf Engl. 2009;25:1754–60.
Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. Twelve years of SAMtools and BCFtools. GigaScience. 2021;10:giab008.
Thai QK, Chung DA, Tran H-D. Canis mtDNA HV1 database: a web-based tool for collecting and surveying Canis mtDNA HV1 haplotype in public database. BMC Genet. 2017;18:60.
Suchodolski JS, Camacho J, Steiner JM. Analysis of bacterial diversity in the canine duodenum, jejunum, ileum, and colon by comparative 16S rRNA gene analysis. FEMS Microbiol Ecol. 2008;66:567–78.
Vital M, Gao J, Rizzo M, Harrison T, Tiedje JM. Diet is a major factor governing the fecal butyrate-producing community structure across Mammalia Aves and Reptilia. ISME J. 2015;9:832–43.
Strauss J, Kaplan GG, Beck PL, Rioux K, Panaccione R, Devinney R, et al. Invasive potential of gut mucosa-derived Fusobacterium nucleatum positively correlates with IBD status of the host. Inflamm Bowel Dis. 2011;17:1971–8.
Ohkusa T, Sato N, Ogihara T, Morita K, Ogawa M, Okayasu I. Fusobacterium varium localized in the colonic mucosa of patients with ulcerative colitis stimulates species-specific antibody. J Gastroenterol Hepatol. 2002;17:849–53.
Bullman S, Pedamallu CS, Sicinska E, Clancy TE, Zhang X, Cai D, et al. Analysis of Fusobacterium persistence and antibiotic response in colorectal cancer. Science. 2017;358:1443–8.
Vázquez-Baeza Y, Hyde ER, Suchodolski JS, Knight R. Dog and human inflammatory bowel disease rely on overlapping yet distinct dysbiosis networks. Nat Microbiol. 2016;1:16177.
Suchodolski JS, Dowd SE, Wilke V, Steiner JM, Jergens AE. 16S rRNA gene pyrosequencing reveals bacterial dysbiosis in the duodenum of dogs with idiopathic inflammatory bowel disease. PloS One. 2012;7:e39333.
Minamoto Y, Minamoto T, Isaiah A, Sattasathuchana P, Buono A, Rangachari VR, et al. Fecal short-chain fatty acid concentrations and dysbiosis in dogs with chronic enteropathy. J Vet Intern Med. 2019;33:1608–18.
Chevrot R, Carlotti A, Sopena V, Marchand P, Rosenfeld E. Megamonas rupellensis sp. nov., an anaerobe isolated from the caecum of a duck. Int J Syst Evol Microbiol. 2008;58 Pt 12:2921–4.
Sakon H, Nagai F, Morotomi M, Tanaka R. Sutterella parvirubra sp. Nov. and Megamonas funiformis sp. nov., isolated from human faeces. Int J Syst Evol Microbiol. 2008;58 Pt 4:970–5.
Lee H, Lee JH, Koh S-J, Park H. Bidirectional relationship between atopic dermatitis and inflammatory bowel disease: a systematic review and meta-analysis. J Am Acad Dermatol. 2020;83:1385–94.
Shi X, Chen Q, Wang F. The bidirectional association between inflammatory bowel disease and atopic dermatitis: a systematic review and meta-analysis. Dermatol Basel Switz. 2020;236:546–53.
Penders J, Gerhold K, Stobberingh EE, Thijs C, Zimmermann K, Lau S, et al. Establishment of the intestinal microbiota and its role for atopic dermatitis in early childhood. J Allergy Clin Immunol. 2013;132:601-607.e8.
Brenner DJ, Fanning GR, Steigerwalt AG, Orskov I, Orskov F. Polynucleotide sequence relatedness among three groups of pathogenic Escherichia coli strains. Infect Immun. 1972;6:308–15.
Devanga Ragupathi NK, Muthuirulandi Sethuvel DP, Inbanathan FY, Veeraraghavan B. Accurate differentiation of Escherichia coli and Shigella serogroups: challenges and strategies. New Microbes New Infect. 2018;21:58–62.
Hong P-Y, Lee BW, Aw M, Shek LPC, Yap GC, Chua KY, et al. Comparative analysis of fecal microbiota in infants with and without eczema. PloS One. 2010;5:e9964.
Penders J, Thijs C, van den Brandt PA, Kummeling I, Snijders B, Stelma F, et al. Gut microbiota composition and development of atopic manifestations in infancy: the KOALA Birth Cohort Study. Gut. 2007;56:661–7.
Lewis S, Nash A, Li Q, Ahn T-H. Comparison of 16S and whole genome dog microbiomes using machine learning. BioData Min. 2021;14:41.
Meason-Smith C, Olivry T, Lawhon SD, Hoffmann AR. Malassezia species dysbiosis in natural and allergen-induced atopic dermatitis in dogs. Med Mycol. 2020;58:756–65.
We thank the dogs and their owners for participating in this study. Authors also thank Ms. Miriam Freitag and Ms. Petra Langenstrassen for their excellent technical assistance. We acknowledge to Dr. Atsushi Yamamoto for supporting communications and for organising clinical data and samples. MT, AK, IW, MO, TL, and HB acknowledge computational support from the OMICS compute cluster at the University of Lübeck. IW and HB acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy-EXC 22167-390884018. We would like to thank the late Dr. Ayako Kamida for supporting clinical samples organisation at the initial phase of the project.
Open Access funding enabled and organized by Projekt DEAL. This work was funded by WALTHAM Foundation grant.
Ethics approval and consent to participate
The Clinical Animal Trials and Research Ethics Committee of the Tokyo University of Agriculture and Technology approved all studies (approval no. 0016021). All participating dog owners were provided with information about the study and gave their written consent before their dogs were included.
Consent for publication
All authors have reviewed the manuscript.
AW is an employee of Royal Canin SAS, Mars Petcare and that the WALTHAM Foundation is funded by Mars Petcare.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. The cohort of Shiba Inu dogs tested in this study (a). Twelve skin swab samples and a faecal sample were collected from each dog in each group, resulted in a total of 12 ´ (10+ 10+ 10 + 20) = 600 skin swab samples and 1 ´ (10+ 10 + 10 + 20) = 50 faecal samples. The twelve skin sites include a buccal sample was collected from three groups; cADpre, cADtreat and Healthy, 1 ´ (10+ 10 + 20) = 40 samples. The twelve skin sites are shown in (b); 1. perilabial area (lips), 2. medial pinnae (concave pinnae), 3. axilla, 4. front paw (dorsal and palmar sides combined), 5. hind paw (dorsal and palmar sides combined), 6. cubital flexor (elbow fold), 7. palmar metacarpal (from carpal to metacarpal pad), 8. flank, 9. inguinal area (groin), 10. abdomen, 11. perineum (from genitalia to anus) and 12. ventral tail (proximal). Age distribution of the participating dogs (c), CADESI-04 (d) and PVAS (e) values for each dog were plotted in the graphs. The values of healthy dogs were given as 0. Asterisks (*) indicate the statistically significant difference; **p < 0.01, ***p < 0.001, ****p < 0.0001 (Kruskal-Wallis test). Ten dogs with untreated cAD dogs (cADpre) received oclacitinib treatment for 2 weeks (cADpost). The treatment efficacy was shown as a statistically significant reduction in CADESI-04 (**p < 0.01, Wilcoxon test; f) and PVAS values (p = 0.0078, Wilcoxon test; g). For these 10 dogs, the CADESI-04 score in 12 different skin sites, the perilabial area, medial pinna, axilla, front paw, hind paw, cubital flexor, palmar metacarpal, flank, inguinal area, abdomen, perineum, and ventral tail, is shown in (h). A significant reduction of CADESI-04 score was observed in inguinal area (*p < 0.05; two-way ANOVA).
Figure S2. For alpha diversity, Shannon indices were compared between four groups: healthy dogs (Healthy), cAD-affected dogs without medical treatment (cADpre), cADpre dogs that received 2 weeks of oclacitinib treatment (cADpost), and cAD-affected dogs already on oclacitinib treatment (cADtreat). The comparisons were made in 12 skin sites: abdomen (a), axilla (b), cubital flexor (c), flank (d), front paw (e), hind paw (f), inguinal area (g), medial pinnae (h), palmar metacarpal (i), perilabial area (j), perineum (k) and ventral tail (l). Asterisks (*) indicate p-values showing significant correlations (p< 0.05) are presented.
Figure S3. Beta diversity was analysed by CAPscale analysis of Aichison distances and was compared between four groups: healthy dogs (Healthy), cAD-affected dogs without medical treatment (cADpre), cADpre dogs that received 2 weeks of oclacitinib treatment (cADpost), and cAD-affected dogs already on oclacitinib treatment (cADtreat) in 12 skin sites: abdomen (a), axilla (b), cubital flexor (c), flank (d), front paw (e), hind paw (f), inguinal area (g), medial pinnae (h), palmar metacarpal (i), perilabial area (j), perineum (k) and ventral tail (l).
Figure S4. Differentially abundant bacterial taxa at the phylum level between two groups; healthy and cADpre, cADtreat and healthy, or cADpre and cADpost in four selected skin sites; abdomen (a), front paw (b), perilabial area (c), and ventral tail (d). Differentially abundant taxa (q < 0.1) are presented in the plots; blue dots refer to phyla identified by corncob, and orange dots refer to phyla identified by corncob and ANCOM-BC.
Figure S5. Differentially abundant bacterial taxa at the genus level between two groups; healthy and cADpre, cADtreat and healthy, or cADpre and cADpost in four selected skin sites; abdomen (a), front paw (b), perilabial area (c), and ventral tail (d). Differentially abundant taxa (q < 0.1) are presented in the plots; blue dots refer to genera identified by corncob, and orange dots refer to genera identified by corncob and ANCOM-BC.
Figure S6. Correlation of pruritis score (PVAS) and bacterial taxa at the phylum level (a) and at the genus level (b) in the four skin sites: abdomen, front paw, perilabial area and ventral tail. P values showing significant correlations (p < 0.05) are presented.
Figure S7. Of the twelve investigated skin sites, the abundance of Staphylococcus was associated with cAD in the skin sites, showing a significant difference at five sites when compared between dogs with untreated cAD (cADpre) and healthy dogs; these sites were the abdomen (a), axilla (b), cubital flexor (c), palmar metacarpal (d), and perilabial area (e). Three sites, the flank (f), front paw (g), and hind paw (h), also showed a tendency towards a higher abundance of Staphylococcus in cADpre dogs than in healthy dogs. Marginal effects for each comparison are shown in the left column plots with 95% confidence intervals (whiskers); the p value for the comparison between cADpre and Healthy is shown. The p values were corrected using the Turkey method for comparing a family of three estimates. The abundance of Staphylococcus in each of these skin sites correlated with its respective CADESI-04 score (i-p).
Figure S8. Data related to Fig. 3. The mean abundances of the 10 predominant Staphylococcus spp. are presented in the 12 skin sites: abdomen, axilla, cubital flexor, flank, front paw, hind paw, inguinal area, medial pinnae, palmar metacarpal, perilabial area, perineum and ventral tail. Abundances are shown per study group; healthy dogs, cAD-affected dogs before treatment (cADpre), cADpre dogs that received 2 weeks of oclacitinib treatment (cADpost) and cAD-affected dogs already on treatment (cADtreat) are presented.
Figure S9. Data related to Fig. 4. Differentially abundant taxa in the gut at the phylum level between the two groups (p < 0.1) are shown. The groups compared are indicated in each panel: (left) cAD-affected dogs that did not receive medical treatment (cADpre) vs. healthy dogs (Healthy); (middle) cAD-affected dogs on medical treatment (cADtreat) vs. Healthy; (right) cADpre dogs that received 2-week oclacitinib treatment (cADpost) vs. cADpre. Blue dots refer to phyla identified by corncob, and orange dots refer to phyla identified by corncob and ANCOM-BC.
Figure S10. Correlation of bacterial taxa in the gut and canine clinical score (CADESI-04; total score; a and b) or pruritis score (PVAS; c and d) at the phylum level (a and c) and at the genus level (b and d). P values showing significant correlations (p< 0.05) are presented.
Figure S11. The PHRED-scaled genotype likelihood of nucleotide variants in the mitochondrial genome identified in each dog are presented in a heatmap. This value ranges between 0 and 255, with 0 being the strongest support of a variant (displayed in black). Positions with no variant call are displayed in white. These variant patterns divide the dogs into three clusters that agree with the mitochondrial haplogroups assigned: A, B and C. The two mitochondrial haplogroups with more than 10 dogs included, haplogroup A (excluding A64, because its variant profile is markedly different) and C, were selected for further analysis of the faecal microbiota.
Figure S12. An association between the mitochondrial haplogroup and skin microbiota in axilla was evaluated. Alpha diversity was compared between skin bacteria in axilla of dogs with haplogroup A and that of dogs with haplogroup C in 3 groups: cADpost, cADtreatment and healthy (a, left panel), and in all three groups together (b, right panel). No difference was observed in alpha diversity between dogs with haplogroup A and those with C in both models (panel a: p = 0.8531; panel b: beta p = 0.8512, estimate = -0.042). Beta diversity was also compared between these two haplogroups (c). The Aitchison distance between haplogroups A and C showed no significant differences (p = 0.7956). Neither haplogroup (p = 0.7971) nor the separation at the x-axis (CAP1, p = 0.7955) was significant; p values were calculated using 99,999 permutations using an ANOVA-like permutation test for constrained correspondence analysis. The differentially abundant bacterial taxa between haplogroups A and C at the phylum level (d) and at the genus level (e) are presented (q < 0.1); blue dots refer to genera identified by corncob, and orange dots refer to genera identified by corncob and ANCOM-BC.
Figure S13. An association between the mitochondrial haplogroup and skin microbiota in medial pinnae was evaluated. Alpha diversity was compared between skin bacteria in medial pinnae of dogs with haplogroup A and that of dogs with haplogroup C in 3 groups: cADpost, cADtreatment and healthy (a, left panel), and in all three groups together (b, right panel). No difference was observed in alpha diversity between dogs with haplogroup A and those with C in both models (panel a, p = 0.7110; panel b, beta: p = 0.6320, estimate = -0.094). Beta diversity was also compared between these two haplogroups (c). The Aitchison distance between haplogroups A and C showed no significant differences (p = 0.3495). Neither haplogroup (p = 0.3497) nor the separation at the x-axis (CAP1, p = 0.3498) was significant; p values were calculated using 99,999 permutations using an ANOVA-like permutation test for constrained correspondence analysis. The differentially abundant bacterial taxa between haplogroups A and C at the phylum level (d) and at the genus level (e) are presented (q < 0.1); blue dots refer to genera identified by corncob, and orange dots refer to genera identified by corncob and ANCOM-BC.
Figure S14. An association between the mitochondrial haplogroup and skin microbiota in front paw was evaluated. Alpha diversity was compared between skin bacteria in front paw of dogs with haplogroup A and that of dogs with haplogroup C in 3 groups: cADpost, cADtreatment and healthy (a), and in all three groups together (b). No difference was observed in alpha diversity between dogs with haplogroup A and those with C in both models (a, p = 0.903; b, beta: p = 0.9210, estimate =-0.028). Beta diversity was also compared between these two haplogroups (c). The Aitchison distance between haplogroups A and C showed no significant differences (p = 0.8515). Neither haplogroup (p = 0.8516) nor the separation at the x-axis (CAP1, p = 0.8481) was significant; p values were calculated using 99,999 permutations using an ANOVA-like permutation test for constrained correspondence analysis. The differentially abundant bacterial taxa between haplogroups A and C at the genus level (d) are presented (q < 0.1); blue dots refer to genera identified by corncob, and orange dots refer to genera identified by corncob and ANCOM-BC. No differentially abundant bacterial taxa at the phylum was identified between the haplogroups.
Table S1. Age and sex distribution of dogs involved in the study. Table S2. Correlation analysis between bacterial taxa and clinical score (CADESI-04) to select the skin sites for further analysis. Table S3. Statistical summary of the differentially abundant skin bacterial taxa at the phylum level. Table S4. Statistical summary of the differentially abundant skin bacteria at the genus level. Table S5. Correlation of Staphylococcus abundance and CADESI-04 score in each skin site was compared between dog groups. Table S6. Mean relative abundance of the 10 most abundant Staphylococcus spp. in the skin samples form healthy, cADpre, cADpost and cADtreat dogs. Data from 12 skin sites are presented. Percent identity of the best hit from NCBI Nucleotide BLAST is presented in each classified species. Table S7. Data related to Fig. 3. Statistical summary of differentially abundant Staphylococcus spp. between dog groups. Table S8. Mean relative abundance of bacterial taxa at the phylum level in faecal samples from all four dog groups. Table S9. Mean relative abundance of bacterial taxa at the genus level in the faecal samples in all four dog groups. Table S10. Statistical summary of the differentially abundant faecal and skin bacterial taxa between mtDNA haplogroups A and C.
About this article
Cite this article
Thomsen, M., Künstner, A., Wohlers, I. et al. A comprehensive analysis of gut and skin microbiota in canine atopic dermatitis in Shiba Inu dogs. Microbiome 11, 232 (2023). https://doi.org/10.1186/s40168-023-01671-2