Low (15 ppb) and high (522 ppb) levels of aflatoxin-contaminated kibble were measured according to standard operating procedures of the Plant Industries Division of the Missouri Department of Agriculture, which uses the AOAC Official method 2005.08: Aflatoxins in Corn, Raw Peanuts, and Peanut Butter (Liquid chromatography with Post-Column Photochemical Derivatization). Samples were subsequently sent to the Center for Veterinary Medicine, FDA for metagenomic analysis. Kibble measuring 15 ppb aflatoxin were considered “low” and kibble measuring 522 ppb were considered “high” level aflatoxin samples for the purposes of this study. Aflatoxin at 15 ppb in kibble is below the FDA action level for aflatoxins of 20 ppb in pet food (20). An additional LCMS Multiple Mycotoxin evaluation was conducted by a third-party Feed Evaluation Lab (Cumberland Valley Analytical Services, Zullinger, PA). Mycotoxins evaluated included aflatoxin, fumonisin, ochratoxin A, and zearalenone in control and contaminated kibble.

“Control” samples were taken from the exact same brand of dog food but different lot number (not part of the recalled aflatoxin-contaminated kibble) to represent a baseline. Kibble ingredients listed chicken by-product meal, corn, wheat, meat meal, rice bran, chicken fat, dried beet pulp, whitefish meal, flaxseed, salt, potassium chloride, choline chloride, vitamins, and minerals. The label described a composition of at least 26% protein, 15% fat with maximum fiber at 6%, moisture at 10%, and 3,645 kcal/kg.

Three 2.5 g portions of each level of aflatoxin-contaminated kibble and control were ground to a fine powder in a Qiagen Tissue Lyser II at 30 Hz (1800 oscillations per minute) for 1.3 min per sample. Replicates of 100 mg of powder from the pooled mixture (7.5 gram) were used for DNA extraction. DNA extraction was performed using the Zymo High Molecular Weight DNA extraction kit according to the manufacturer’s specifications including an extra lysing step using PBS and lysozyme (21). Libraries of DNA were created using the Illumina DNA Library Prep Kit according to the manufacturer’s specified protocols and sequencing was performed on a NextSeq 2000 using a high-throughput kit according to previously described methods (22). An additional NextSeq 2000 high-throughput sequencing run was performed with only two replicates of low- and high-level samples to achieve a sequencing depth between 100 and 250 million reads per replicate to evaluate how increased read depth impacted incidence of key species.

Sequence data were screened for quality metrics using Trimmomatic (23) and analyzed using in-house FDA pipelines and databases as previously described (22). Determination of bacterial and fungal composition from shotgun sequencing was conducted using custom C++ programs developed to compile a k-mer signature database containing multiple unique 30 bp sequences per species and then identify each read in the input file using the 30 bp probes. For each bacterial or fungal species or subspecies, each non-duplicated 30-mer from a reference whole genome sequence was placed into a database. Any k-mers not found in at least 2/3 of a set of additional genome sequences of the same species were removed and k-mers found in genomes of other species were removed. Normalization was performed to correct for bias due to differing number of k-mers per database entry and results were tabulated as percent of identified reads (contribution to the microbial population of identified species) for each database entry. Results can be expressed by raw hits or as relative abundance. Both are useful when detecting low abundance organisms in complex metagenomes. Annotation was also done using the Cosmos ID cloud-based application (CosmosID Metagenomics Cloud, app.cosmosid.com, CosmosID Inc., www.cosmosid.com) with Fungal Database Version 1.2 (24). LEfSe (Linear discriminant analysis effect size) by the Huttenhower biobakery¹ (25) was calculated in the COSMOSID application. LEfSe (Linear discriminant analysis effect size) is an algorithm used for biological biomarker discovery. Features such as genes, pathways, or taxa were identified for each treatment and the non-parametric factorial Kruskal-Wallis (KW) sum-rank test (26) was used to identify significant differential abundance of specific features between treatments. Linear Discriminant Analysis (LDA) was used to estimate the effect size of each differentially abundant feature and rank the feature accordingly. For the bar chart presented in the results section, the LDA ranged from 2.9 to 5.27 and the p value ranged from 0.022 to 0.035. Sequence data annotations were visualized using graphs created by the R Tidyverse package. ²

Metagenomic reads were mapped with Diamond (v2.0.5) (27) BLASTX (≥95% identity and those ≥90% read coverage) to a database of genes involved in mycotoxin biosynthetic pathways: aflatoxin, deoxynivalenol, nivalenol, ochratoxin, patulin, sterigmatocystin, T2_toxin, and tenuazonic_acid. Genes were identified and downloaded from the MetaCyc database (28). ³ Identity of reads aligning to mycotoxin genes was confirmed by BLASTX, aligning them to the NCBI nr database online.

To describe composition and relative abundance of ingredients, metagenomic data were created for three replicates of control, three replicates of low-concentration (15 ppb), and three replicates of high-concentration (522 ppb) aflatoxin-contaminated kibble. An average of 182 million sequences per replicate were used in downstream analyses. The relative abundance of macro ingredients [annotated by mitochondrial DNA (29) and occurring at greater than 1% of normalized data] included Zea (corn), Gallus (chicken), Triticum (wheat), Soya, (soybean), Bos (cow, ox, bull, yak, cattle), and yeast. Further refinement of species annotations described Z. mays, G. gallus, G. gallus ssp. spadiceus, T. aestivum, B. taurus, several species of Saccharomyces (cerevisiae, pastorianus, and pastorianus Weihenstephan), and Fusarium verticillioides.

Evaluation of the metagenomic data from control, low-, and high-level aflatoxin-containing kibble correlated with relative abundance of DNA of Aspergillus species, i.e., the low-level aflatoxin-contaminated kibble had a low abundance of DNA from Aspergillus species and the high-level kibble had a higher relative abundance of Aspergillus species. The incidence of Aspergillus or any other taxon can be expressed by raw hits to the genome of interest or as relative abundance after normalization by k-mer count per organism and number of sequence reads. Almost no Aspergillus DNA was detected in controls (Figure 1). The Aspergillus species identified in control kibble was predominantly A. glaucus. Figure 1A shows the top twenty most abundant Aspergillus species identified in control, low- (15 ppb), and high- (522 ppb) levels of aflatoxin-contaminated dog kibble using relative abundance normalization. Figure 1B shows the raw k-mer hits to the top 20 Aspergillus species seen in the metagenomes. Both approaches, relative abundance and raw k-mer hits, are useful when attempting to detect low abundance DNAs from low abundance organisms. There was an extensive taxonomic range of Aspergillus species identified in the kibble. Production of aflatoxin may have been associated with a single species or potentially multiple species. Aspergillus oryzae, flavus, and phoenicis were the most abundant species observed in high-level (522 ppb) aflatoxin-contaminated kibble. Aspergillus flavus and A. parasiticus are two of the most traditionally recognized producers of aflatoxin in pre- and post-harvest commodities (12). These species likely played a role in aflatoxin production in the kibble as A. oryzae and A. phoenicis are not known to produce aflatoxins (30, 31).

While aflatoxin B1 was the primary focus of the chemical evaluation of the kibble due to its acute toxicity, concentrations of fumonisins (B1,B2,B3), and ochratoxin A were also detected by LCMS (Table 1).

Fusarium verticillioides is reported to produce fumonisins and zearalenone and was observed in both low and high samples of kibble but not in controls (Figure 2A). Additional Fusarium species, including F. annulatum, proliferatum, dlaminii, siculi, irregulare, pilosicola, and globosum, were primarily associated with contaminated kibble and were not observed in controls (Figure 2A). Deoxynivalenol is also produced by Fusarium species but was not measured in this study.

Penicillium spp. Are reported to produce Ochratoxin A and patulin (32). Penicillium oxalicum, subrubescens, brasilianum, ochrochloron, roqueforti, and citrinum were primarily observed in contaminated samples whereas P. verrucosum, nordicum, and freii were the most abundant species observed in controls (Figure 2B). It is possible that P. oxalicum or P. subrescens played a role in the observed ochratoxin levels (23 ppb) in ‘high’ samples. However, P. oxalicum was also observed in low-level samples for which no ochratoxin was quantified, so perhaps other genera that were differentially enriched in high-level samples were responsible for the ochratoxin observed by LCMS (Table 1).

Linear discriminant analysis of fungal species in low- and high-level contaminated kibble was calculated using the Huttenhowever biobakery tools (25) available in the COSMOSID analysis pipeline (CosmosID Metagenomics Cloud, app.cosmosid.com, CosmosID Inc., www.cosmosid.com). This tool functions as a biomarker discovery tool, comparing differential abundance of taxonomy or functional genes in different treatment groups. Used with fungal taxonomy for control, low-, and high-level aflatoxin contaminated kibble, it was clear that high-level kibble contained a greater abundance of Aspergillus species than low-level kibble (Figure 3). Interestingly, control kibble had significantly more Saccharomyces species than the high-level kibble (p < 0.022 to 0.035). It is common practice to add yeast to ingredients that contain aflatoxin. Yeast species have been shown to bind aflatoxin B and, thus, reduce the toxic impact of consumption of aflatoxin-contaminated food (33). This approach has been shown to provide a protective effect to broiler chickens consuming aflatoxin B1-contaminated food (34). Yeast was listed as an ingredient in all dog kibble.

Additionally, a greater relative abundance of Alternaria species were observed in low-level aflatoxin-contaminated kibble. While Alternaria species were observed in metagenomic data, Alternaria-associated toxins (altenuene, alternariol, alternariol monomethylether, tentoxin, and tenuazonic acid) were not measured in this study. Figure 4 provides a breakdown of the species of Alternaria across each kibble type.

For Alternaria, species composition from control to contaminated kibble was more stable than what was observed for genera like Aspergillus, Penicillium, and Fusarium (Figures 1, 2). Fingerprints of differential abundance of fungal taxa (Figure 5) could potentially be used to quickly detect lots of feed or ingredients that may be at higher risk of aflatoxins or other mycotoxin exposure. When Saccharomyces species (used to bind aflatoxin) have been depleted, there is likely a greater risk of toxicity exposure.

Beyond identification of the taxonomic structure of kibble, metagenomic data also described functional genes involved in mycotoxin production and thus provided further confirmation of mycotoxin pathways such as deoxynivalenol, nivalenol, and aflatoxin. Not surprisingly, the most numerous hits to genes involved in the aflatoxin production pathway were seen in the high-level aflatoxin-contaminated samples (Table 2). The doublet and singleton hits to nivalenol (NIV) and deoxynivalenol (DON) in control samples may be indicative of low levels of NIV or DON in these samples.

Stenocarpella maydis, another potential toxin producer, was observed in the contaminated kibble samples but not in controls. Stenocarpella was identified as a biomarker for low-level aflatoxin samples where it was significantly enriched (Figure 3). Toxic metabolites including diplodiatoxin, chaetoglobosins K and L, and (all-E)-trideca-4,6,10,12-tetraene-2,8-diol associated with Stenocarpella are not commonly evaluated and simple ELISA or other HPLC based tests are not readily available.

Bacterial profiles for controls and aflatoxin-contaminated dog kibble are shown in Figure 6. Three replicates of control, low-, and high-level aflatoxin-contaminated dog kibble were averaged and taxa occurring at greater than 3% of total data were graphed to summarize bacterial features of each kibble. Enterobacter, Serratia, and Kosakonia genera were observed in toxin-contaminated kibble and not in control kibble.

Control kibble had different antimicrobial resistance genes than those observed in aflatoxin-contaminated kibble. This is most likely driven by the different bacteria that were present in contaminated samples contrasted to controls. Genes involved in macrolide resistance were seen in contaminated kibble as well as beta-lactam genes bla R1 and bla I, which influence expression of blaZ and penicillin binding protein 2a (PBP 2a) and contribute to methicillin resistance (35, 36). The bla 1 gene was observed in control samples, typically associated with ampicillin resistance.