Fecal samples from two BFFs exhibiting symptoms of coccidiosis were collected over two consecutive days (“Ferret1” and “Ferret2”). Microscopic slides of fecal samples were prepared using both a direct fecal exam and Sodium nitrate solution (Fecasol^®, Vetoquinol USA, Inc. Fort Worth, TX, USA). Coccidia oocysts were visualized under a compound microscope (Abaxis, 3000-LED Series) at 400 magnification (40x objective, 10x eyepiece). Fecal samples were collected in DNA preservative buffer (DNA/RNA Shield TM, Zymo Research Corp. Cat. No. R1108, Irvine, CA, USA) and submitted to MiDOG LLC for DNA extraction and sequencing. Sequencing libraries were prepared with 10bp unique dual indexes and an Illumina DNA Prep Kit (Illumina, San Diego, CA, USA) following the manufacturer’s protocol. Each library was quantified using Qubit (Thermo Fisher Scientific) and pooled together at equal abundance, and the final pool was quantified using qPCR. The metagenomic DNA was sequenced using the Illumina Novaseq X platform, generating 151bp paired end reads. Raw sequence reads were trimmed to remove-low quality fractions and adapters using Trimmomatic v.0.33 [15]. Low-quality fractions were removed using a sliding window with a 6bp window size and a quality cutoff of 20. Adaptors and reads smaller than 70bp were also removed. After that, host-derived reads were removed using Kraken2 [16] against common eukaryotic host genomes. Sdust [17] was used to detect and remove the low-diversity reads.

Microbial profiles for each sample were generated using Sourmash v.4.8.11 [18] with the sketch command, using a k-mer size of 51 and a scaling factor of 1000. The resulting 51-bp k-mers were compared against a collection of curated and decontaminated genomes from bacteria, fungi, and protists to identify the microbial taxa present in the samples. Reference databases included the GTDB species representative database (RS207), GenBank pre-formatted databases (v. 2022.03), and internally curated MiDOG reference database. Together, these databases encompass 64,829 bacteria, 10,291 fungal,1,189 protozoan, 879 metazoans (nematodes and platyhelminths), and 47,286 viral genomes. Only species with total kmer coverage exceeding 5 kilobases (kb) were retained to remove potential false positives. Reads were mapped back to the genomes identified by Sourmash using BWA-MEM [19], and microbial abundance was determined based on the counts of mapped reads. Sample signatures were compared using the Sourmash compare command employing angular similarity as the similarity metric. The Sourmash gather function was subsequently used to estimate the proportion of each query signature that match the reference databases. Taxonomic annotation of the gathered results was performed using the Sourmash tax annotate. The resulting outputs from the sourmash analyses were imported into R for downstream visualization and statistical analysis. Bar plots showing the proportions of query matches across different taxonomic groups and upset plots indicating the microbial taxonomic lineages present in each sample were rendered using the sourmashconsumr package [20]. Hierarchical visualization of taxonomic relationships was conducted using the metacoder package [21]. Taxonomic annotation data were converted to a MetaCoder object via the sourmashconsumr::from_taxonomy_annotate_to_metacoder function. The Compare_groups function was used to compare infected and uninfected groups. Outputs from sourmash_tax_annotate were converted into phyloseq objects using the sourmashconsumr package. Rarefaction curves were generated using phyloseq v.1.52.0 [22] and Vegan v.2.7-1 [23] packages in R, with the parameter step = 10000.

Alpha diversity was assessed using the “Observed” and “Shannon” diversity indices as implemented in the phyloseq package. Statistical significance in alpha diversity between groups was evaluated using the Wilcoxon rank-sum test, with a p-value < 0.05 considered statistically significant. Beta diversity was quantified using Bray-Curtis dissimilarity and visualized through a principal coordinate analysis (PCoA) plot generated with the vegan package. Differences in community composition were tested using permutational multivariate analysis of variance (PERMANOVA) (ADONIS) with 999 permutations in the vegan package.

Differentially abundant microbial species in Coccidia-positive and -negative fecal samples were identified using the DESeq2 package [24]. Pairwise comparisons were performed between the Coccidia “Positive” and “Negative” groups, with the significance threshold (α) set to 0.01. Log₂ fold changes in abundance were calculated to quantify differential enrichment, and the results were visualized using the ggplot2 package in R.

Fecal samples of two BFFs diagnosed with coccidiosis were collected over two consecutive days. Coccidia oocysts were only detected by microscopy on Day 2 (Fig. 1A). All fecal samples collected were subjected to metagenomic sequencing. The number of raw reads generated for each sample is summarized in Table I. Sourmash with a 51bp k-mer size and a scaling factor of 1000 was used for compositional analysis and taxonomic profiling of the microbial community. The 51bp k-mers were queried against the GTDB species representative database (RS207), pre-formatted GenBank database (v. 2022.03), and curated MiDOG database to assign taxonomic classifications. Samples collected on the first day (Coccidia-negative) contained 75,158 and 2,550,412 classified reads, whereas samples collected on the second day (Coccidia-positive) contained 1,782,890 and 418,045 classified reads, respectively (Table I). Overall, the Coccidia-negative and -positive groups exhibited comparable mean classified read counts (Table II).

Across all samples, 51bp kmers were searched against three databases to determine taxonomic assignments. Among the classified reads, the majority were identified as bacterial, comprising 49%, 74%, 76%, and 39% of each respective sample (Fig. 1B). Each sample also contained at least 25% unclassified reads, suggesting the presence of novel or uncharacterized microbial species that were not represented in these databases.

Rarefaction curves were generated for each sample by subsampling 10,000 reads per iteration (Supplementary Fig. 1). The curves for Ferret2_Coccidia_Neg, Ferret1_Coccidia_Pos, and Ferret2_Coccidia_Pos samples plateaued, indicating that the sequencing depth was sufficient to capture most of the species present. In contrast, the curve for the Ferret1_Coccidia_Neg sample did not plateau and terminated near the beginning of the plot, suggesting inadequate sequencing depth to fully represent the microbial diversity of that sample.

A total of 15 unique lineages at the species level were detected across all four samples, seven belonging to the phylum Proteobacteria and four belonging to Actinobacteriota (Fig. 1C). The Ferret2_Coccidia_Neg sample contained 32 species unique to the sample, 18 of which belonged to the phylum Proteobacteria. Notably, six species were shared between the two Coccidia-positive samples, while no species were shared between the two Coccidia-negative samples.

E. coli accounted for approximately 60% and 80% of the microbial abundance in the two Coccidia-negative samples, respectively (Supplementary Fig. 2). Upon Coccidia shedding, the relative abundance of E. coli in Ferret 1 decreased from 58% to 28%, while Bacterioides fragilis remained stable at approximately 30%. Conversely, Fusobacterium sp. and Enterococcus faecalis increased to approximately 10%. For Ferret 2, the relative abundance of E. coli decreased from 80% to 50%, whereas that of both Fusobacterium sp. and Peptostreptococcus russellii increased to approximately 17%.

At the phylum level, Proteobacteria had the highest number of lineages across the samples (Supplementary Fig. 3). The family Enterobacteriaceae and the genera Escherichia and Enterococcus exhibited the greatest lineage diversity. Members of Firmicutes_A, including Clostridia, Peptostreptococcaceae, Bacteroidales, and Fusobacterium, were significantly enriched in the coccidia-positive group. In contrast, Enterobacteriaceae, including Escherichia and Hafnia were more abundant in the coccidia-negative group. Yersinia pestis, the causative agent of sylvatic plague, appeared was enriched in the Coccidia-positive group; however, its relative abundance compared to the overall microbiome was low, 0.0004% in Ferret 2_Coccidia_Neg and 0.0193% in Ferret 1_Coccidia_Pos.

Only one protozoan species was detected exclusively in the two Coccidia-positive samples, corroborating microscopy-based diagnosis. Taxonomic assignment identified this protist as belonging to the genus Cystoisospora, which includes a known pathogen responsible for coccidiosis in BFFs A total of 105 classified bacterial species and one protozoan species were identified across all samples. Both groups exhibited similar median observed species richness for total and bacterial taxa, with no statistically significant differences (p-value = 1.0 and 1.0, respectively) (Fig. 2A and 2C). Although the Shannon diversity index was slightly higher in the Coccidia-positive group for both total and bacterial taxa, these differences were not statistically significant (p-value = 0.33 and 0.33, respectively) (Fig. 2B and 2D). Despite detection of the protozoan species only in the Coccidia-positive group, there were no significant differences in the observed (p-value = 0.1) or Shannon diversity indices (p-value = NA) between the Coccidia-positive and -negative groups (Fig. 2E and 2F). PCoA of Bray-Curtis dissimilarity yielded an F-score of 0.86 (p-value = 1.0) for total beta diversity and 0.95 (p-value = 0.67) for bacterial beta diversity (Supplementary Fig. 4), indicating that between-group dissimilarity was not significantly greater than within-group dissimilarity.

The relative abundance of microbial taxa was compared across the Coccidia-positive and -negative groups (Fig. 3). In the Coccidia-negative group, members of the phylum Proteobacteria accounted for 91.5% of the total abundance. In contrast, in the Coccidia-positive group, Proteobacteria decreased to approximately 39%, Bacteroidota increased to 26.8%, and both Firmicutes (11.4%) and Fusobacteriota (14.1%) increased to over 10% (Fig. 3A). At the genus level, Escherichia represented 86% and Hafnia represented approximately 8.2% of the total abundance in the coccidia-negative group. In the Coccidia-positive group, Escherichia decreased to approximately 39%, Bacteroides increased to 27%, Enterococcus to 12%, and Fusobacterium to 14.5%, each exceeding 10% relative abundance (Fig. 3B). Among the protozoa, Cystoisospora sp. was the only species detected, appearing exclusively in the Coccidia-positive group. This organism corresponds to Coccidia. which is responsible for the infection. No protozoan species were detected in the coccidia-negative group (Fig. 3C).

DESeq2 analysis revealed that only one bacterial species, Peptostreptococcus russellii (phylum Firmicutes), was significantly enriched in the Coccidia-positive group based on log₂ fold change, with an adjusted p-value cutoff of 0.01 (log2FoldChange = 17.73, adj. p-value = 1.13e-04) (Fig. 4). Conversely, several taxa were significantly depleted in the day-two samples, including Cutibacterium acnes (log2FoldChange = -24.12, adj. p-value = 4.65e-07) and six species from the phylum Proteobacteria, namely Cedecea michiganensis (log2FoldChange = -25.23, adj. p-value = 3.62e-06), Klebsiella planticola (log2FoldChange = -23.52, adj. p-value = 1.38e-05), Escherichia coli (log2FoldChange = -26.15, adj. p-value = 2.15e-06), Escherichia marmotae (log2FoldChange = -25.93, adj. p-value = 2.15e-06), Hafnia alvei (log2FoldChange = -30.00, adj. p-value = 3.90e-08), and Hafnia proteus (log2FoldChange = -24.83, adj. p-value = 4.56e-06).