A comprehensive analysis of gut and skin microbiota in canine atopic dermatitis in Shiba Inu dogs

Background 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. Results 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. Conclusions 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. Video Abstract Supplementary Information The online version contains supplementary material available at 10.1186/s40168-023-01671-2.


Background
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 [7].cAD has a strong association with breed, although the breed predisposition varies among geographical locations [8], 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 [9], and changes in microbial composition, namely, dysbiosis, are known to drive the pathology of complex diseases [10].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 [15].Furthermore, a number of studies have demonstrated that gut dysbiosis is accompanied in human AD patients (reviewed in [16]).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 [17].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 [17]).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 [20].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 [21] and one breeding colony of Shiba Inu dogs [22], or a study of multiple dog breeds [23], 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 [29].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 [30].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) [31].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.

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-NGS PREP 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).

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 [34], with small modifications.In brief, two primer sets [CLF_mtDNA_1_F (TTC TTC GGC GCA TTC CAC AA) and CLF_mtDNA_1_R (GGC ATG CCT GTT AAT GCG AG); CLF_mtDNA_2_F (TCA CTT TAT GCT TAG GGG CCA) and CLF_mtDNA_2_R (CAA CAT TTT CGG GGT ATG GGC)] 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.).

Alpha and beta diversity
In faecal samples, alpha diversity was estimated using DivNet [43] (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 [44] 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 [45]).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 [46] 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 (H 0 ~ sex, H 1 ~ group + sex).For the haplogroup analysis, disease group was incorporated into the NULL model to correct for these effects (H 0 ~ group, H 1 ~ group + haplogroup).To verify the findings, ANCOM-BC [47] (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.
Balances with CADESI-4 or PVAS as the target variable were calculated using a forward-selection method as implemented in selbal [48] (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/ busch lab/ 2022_ canine_ atopic_ derma titis_ paper.

Mitochondrial genome analysis
Raw sequencing data were mapped with bwa (version 0.7.17) [49] against dog reference sequence U96639.2(https:// www.ncbi.nlm.nih.gov/ nucco re/ U96639), and variant calling was performed with bcftools (version 1.15) [50] 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.vnbio logy.com) [51] for haplogroup assignment.The MT analyses workflow used is available at https:// github.com/ TLenf ers/ multi speci es_ mitoc hondr ial_ varia nt_ analy sis.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/ iwohl ers/ 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 Fig. 1 Taxonomic classification of 12 skin sites from healthy dogs, dogs with untreated cAD (cADpre), cAD-affected dogs that received 2 weeks of oclacitinib treatment (cADpost) and cAD-affected dogs that had already been on the treatment (cADtreat) at the phylum level (a) and the genus level (b).The abundance of genus Staphylococcus in 12 skin sites in all 4 dog groups was presented in (c).*p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001.Two-way ANOVA 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 crossvalidation 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 (R 2 > 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 cADaffected 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 cADaffected 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).3 The abundances of the 10 predominant Staphylococcus spp.were compared between the two groups: cADpre vs. healthy, cADtreat vs. healthy and cADpost vs. cADpre at the 12 skin skites: 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).The beta coefficient represents the coefficients obtained from the ANCOM-BC log-linear (natural log) model, and the standard error of the coefficients are shown as whiskers.The x-axis indicates log2 fold-change Fig. 4 Faecal microbiota analysis.Mean relative taxonomic abundances at the phylum level (a) and at the genus level (b) in each dog group are presented.The differentially abundant bacterial genera (q < 0.1) are presented (c), where blue dots refer to genera identified by corncob and orange dots to genera identified by corncob and ANCOM-BC.The comparison of the abundance was conducted between the two groups 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 dogs; right, cADpre dogs that received 2 weeks of oclacitinib treatment (cADpost) vs. cADpre.Alpha diversity was analysed using the DivNet estimate of Shannon, as samplewise estimates (d) and communitywise estimate (e).Beta diversity was evaluated using constrained analysis of principal coordinates of the Aichison distance, demonstrating clear differences between the four groups (adonis: p < 0.001) (See figure on next page.)Fig. 4 (See legend on previous page.)

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).
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.

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 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 [31], 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 (r 2 ) 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.vnbio logy.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).Fig. 5 An association between the mitochondrial haplogroup and gut microbiota was evaluated.Alpha diversity was compared between faecal bacteria 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).When comparing haplogroups A and C in all three groups (regardless of health status), alpha diversity in dogs with haplogroup A showed a higher trend than that in dogs with haplogroup C (p = 0.074).Beta diversity was also compared between these two haplogroups (c).The Aitchison distance between haplogroups A and C demonstrated a difference with p = 0.078.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 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.

Discussion
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 [52].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 [53].In humans, Fusobacterium is known to be highly abundant in various gastrointestinal diseases, such as inflammatory bowel disease (IBD) [54], ulcerative colitis [55], and colorectal cancer [56], suggesting that the impact of Fusobacterium on health and disease in dogs seems to be different from that in humans [57].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 [23].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 [21].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 [58].In dogs, the abundance of Fusobacterium is positively correlated with faecal concentrations of butyrate and propionate [59], 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 [64] 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 [66].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 cADaffected 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 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 [69].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 [70].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.

Conclusions
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.
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 ANOVAlike 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.
Additional file 15: 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.

Fig. 2
Fig. 2 Spearman's correlation of clinical score (CADESI-04) and bacterial taxa at the phylum level (a) and at the genus level (b) at the four skin sites: abdomen, front paw, perilabial area, and ventral tail.P values showing significant correlation (p < 0.05) are shown

Table 1
Statistical summary of the differentially abundant faecal bacterial taxa at the phylum level

Table 2
Statistical summary of the differentially abundant faecal bacterial taxa at the genus level