Enriched Shotgun Sequencing to Assess the Effects of Interventions to Reduce Antimicrobial Use in Neonatal Dairy Calves

Abstract

ABSTRACT

Background

Neonatal diarrhea accounts for 20%–25% of morbidity among calves, and antimicrobial drugs (AMDs) are often administered for treatment. Systematic approaches that mitigate antimicrobial use (AMU) can be effective in decreasing antimicrobial resistance (AMR).

Hypothesis/Objectives

To determine the effects of an algorithmic farm‐based intervention that reduced AMU for diarrhea on the community structure of antimicrobial resistance genes (ARGs) identified in the feces of healthy dairy calves.

Animals

Thirty‐one fecal dairy calf samples collected at two timepoints and farms (N = 7–8 per sampling point) were used. Samples were obtained before AMU reductions and 12 months afterward.

Methods

Target‐enriched shotgun sequencing was performed to characterize all ARGs in samples. Bioinformatics processing and statistical analysis were performed using the AMR++ pipeline, MEGARes AMR database, and R.

Results

Pre‐intervention comparisons showed increased relative abundances (RA) consistent with the AMU on each farm. Intra‐farm results showed that on Farm 1, there were significant increases in the RA of ARGs for tetracyclines (22.1%–27.4%, q = 0.02) and fluoroquinolones (0%–0.1%, q < 0.0001) in the Post period. On Farm 2, significant decreases were seen over time in the RA of ARGs for sulfonamides (9.6%–5.1%, q = 0.006) and fluoroquinolones (0.77%–0.12%, q = 0.004).

Conclusions and Clinical Importance

Despite similar reductions in AMU on both farms, implementing an antimicrobial stewardship algorithm was associated with differing effects on and changes to the fecal resistome.

Abbreviations

• AMD(s)
• antimicrobial drug(s)
• AMR
• antimicrobial resistance
• AMU
• antimicrobial use
• ANOSIM
• analysis of similarities
• ARG(s)
• antimicrobial resistance gene(s)
• CAP
• cationic antimicrobial peptides
• CSS
• cumulative sum scaling
• MDR
• multi-drug resistance
• MLS
• macrolides, lincosamides, and streptogramins
• NMDS
• non-metric multidimensional scaling
• RA(s)
• relative abundance(s)

Introduction

Neonatal calf diarrhea accounts for 20%–25% of morbidity and 50% of mortality in dairy calves < 30 days of age. Antimicrobial drugs (AMDs) are commonly administered with the goal of decreasing adverse health outcomes and promoting animal welfare [1–6]. The fecal microbiota of dairy cattle is well documented [7–11], and the effects of AMU on the fecal/gut microbiota of dairy calves have been similarly noted [12–15]. Implementation of approaches consisting of farmer education in the assessment of calf health and modification of management practices for disease prevention can significantly reduce AMU without reducing farm productivity [4].

Several studies have investigated the relationship between AMU and the recovery of antimicrobial resistance genes (ARGs) using indicator species (Escherichia coli and Salmonella) as a proxy for the broader resistome, with contradictory results [16–19]. Two studies have proposed modifications to AMU policy on-farm as a method for reducing antimicrobial resistance (AMR) in dairy cattle, including limiting the classes of AMDs available for treatment and eliminating combination therapies as well as restricting AMU of critically important AMDs [20, 21]. In both instances, only E. coli was used for evaluation of ARGs. However, microbial populations are complex, with many interactions among bacterial species. The use of indicator organisms precludes evaluation of ARGs in the broader microbial population context and might not give a complete understanding of AMR prevalence due to their single-source nature. Conversely, shotgun metagenomic sequencing allows for more comprehensive coverage of ARGs and genomic data, including from difficult-to-cultivate organisms. Despite these advantages, the use of shotgun metagenomics for the characterization of the fecal resistome of dairy calves is limited [22–24]. Target-enrichment of ARGs expands on the utility of shotgun metagenomics in the evaluation of the resistome as it selectively amplifies the < 1% of ARG-related DNA, allowing for the characterization of even rare resistome elements [25]. To date, no study has used target-enriched metagenomics to evaluate AMU policy changes on the fecal resistome of dairy calves.

Significant reductions (average of 75%) in AMU for the treatment of neonatal calf diarrhea after implementation of on-farm interventions are possible without adverse health outcomes [15]. The Canadian Veterinary Medical Association has since adopted this treatment protocol as part of their prescription guidance for dairy calf diarrhea, and with modifications for beef calves [26, 27]. As such, the outcomes of this policy change on the fecal resistome of these same farms [15] should be examined. The objectives of this study were to use target-enriched shotgun metagenomics to investigate the effects of this farm-based intervention [15] for reduced AMU on the community structure of the fecal resistome of healthy neonatal dairy calves on two farms.

Materials and Methods

Ethics Statement

The samples in this study were obtained as part of a previously published study (University of Guelph Animal Care Committee: eAUP 3793) which assessed the effect of an on-farm decision-making algorithm to reduce AMU for neonatal calf diarrhea [15].

Study Samples

Fecal samples were collected per rectum from clinically healthy dairy calves (1–30 days old) housed on two large commercial dairy farms in south-central Ontario in 2015 and 2016 as part of a study evaluating antimicrobial use in calves after an antimicrobial stewardship intervention [15]. Both farms consisted of Holstein calves that received calf starter rations, oral bovine coronavirus, and E. coli K99 antibodies, and cows that were vaccinated against bovine rotavirus and bovine coronavirus 4 weeks before calving. Some additional management information is listed in Table 1. Further management variables for each of the farms were previously reported [15].

TABLE 1 Farm and management variables on two dairy farms before (Pre) and after (Post) implementation of an AMD treatment protocol.

	Farm 1 Pre	Farm 1 Post	Farm 2 Pre	Farm 2 Post
Calves enrolled (N)	7	8	8	8
Calves from external sources?	Yes	Yes	No	No
Housing	Group	Group	Individual	Individual
Bedding	Sawdust	Sawdust	Shavings	Shavings
Colostrum	4 L in first 4 h	5 L in first 4 h	6 L in first 6 h	6 L in first 6 h
Diet	Pasteurized milk (non-antibiotic)	Milk replacer (non-medicated)	Milk replacer (non-medicated)	Pasteurized milk (non-antibiotic)
Volume/feeding (% of body weight)	15%	15%	12%	12%
Feeding method	Robot machine	Robot machine	Bucket	Bucket
Isolation of sick calves?	No	Yes	No	Yes
NSAIDs	Yes	Yes	Yes	Yes
AMD treatment protocol (diarrhea)	Spectinomycin 30 mg/kg IM every 24 h for 10 days + Lincomycin 15 mg/kg IM every 24 h for 10 days + TMS1 16 mg/kg IM every 24 h for 5 days	TMS1 16 mg/kg IM every 24 h for 3 days	TMS2 1920 mg PO once + Ceftiofur sodium 2.2 mg/kg SC every 24 h for 3 days or TMS1 16 mg/kg IM every 24 h for 3 days	TMS1 16 mg/kg IM every 24 h for 3 days or Ceftiofur sodium 2.2 mg/kg SC every 24 h for 3 days
Calves treated with AMDs (diarrhea)	95%	20%	87%	15%
AMD treatment options for pneumonia (1st line—4th line)	Tulathromycin Florfenicol Danofloxacin Enrofloxacin	Tulathromycin Florfenicol Danofloxacin Enrofloxacin	Florfenicol Tulathromycin CCFA Enrofloxacin	Florfenicol Tulathromycin CCFA Enrofloxacin

Fecal samples were collected from clinically healthy calves (1–30 days of age) on these farms 1–6 weeks before (Pre) implementation of an intervention consisting of the education of farmers on calf health evaluation, refining disease prevention strategies, and the implementation of an algorithm for AMD use in diarrheic calves [15]. Briefly, AMDs were administered for the treatment of diarrhea if calves had a fever (rectal temperature > 39.2°C) or hematochezia alongside changes in demeanor and milk intake. Samples were then collected from a different cohort of clinically healthy calves (1–30 days of age) 12 months after implementation of the algorithm (Post). Calves were excluded from sampling if they had an earlier episode of diarrhea, had other concurrent diseases, or had received any antimicrobial treatment. On both farms, calves that developed diarrhea within 10 days of sampling were excluded, and new calves enrolled due to the potential of microbial alteration associated with prodromic signs of diarrhea. The fecal samples were stored at −80°C from the time of collection until processing. This study's samples were a subsample of the total samples collected for the study, which assessed the effect of an on-farm decision-making algorithm to reduce AMU for neonatal calf diarrhea [15], with DNA yield being the main factor in their inclusion.

Both farms were large commercial dairy farms with free stall housing and automated milking systems. Herd size was approximately 600 cows on Farm 1, with an average milk production of 10 300 kg/cow/year, while Farm 2 had a herd size of approximately 700 cows and had an average milk production of 10 400 kg/cow/year. While calves treated for pneumonia were excluded from this study, it should be noted that several AMDs were used on each farm for the treatment of respiratory diseases (Table 1).

DNA Extraction, Bait-Capture and Enrichment Library Preparation, and DNA Sequencing

As previously described, fecal samples were processed for total DNA extraction in 2015 and 2016, using the standard protocol from the ENZA Stool DNA kit (Omega Bio-Tek, Doraville, GA) [15]. DNA yield was measured with a Qubit 2.0 fluorometer (Invitrogen, Waltham, MA, USA) using a High-Sensitivity assay (1 μL of DNA in 200 μL).

SureSelect^{XT HS} Target Enrichment System for Illumina Paired-End Multiplexed Sequencing Libraries with a custom 0.5–2.9 Mb kit (Agilent Technologies, Santa Clara, CA, USA) was used; with some modifications (Appendix S1) for library preparation to improve on-target sequencing yield for ARGs. A customized bait design targeting 3824 ARGs present in the MEGARes v. 1.0 database [25, 28] was applied. DNA sequencing was performed at the University of Colorado Genomics Shared Resource Facility, using a NovaSeq 6000 (Illumina, San Diego, CA, USA) with a 2 × 150 base-pair paired-end chemistry.

Bioinformatics and Statistical Methods

Sequencing data were processed using the AMR++ v. 1.0 bioinformatics pipeline in conjunction with the MEGARes v. 1.0 resistance gene database [28, 29]. Briefly, read trimming and quality filtering was executed using Trimmomatic [30]. Reads matching the Bos taurus [31] genome were removed using the Burrows-Wheeler Aligner (BWA) and SamTools [32–34]. Non-host reads were subsequently aligned to the MEGARes v. 1.0 database with BWA, and read deduplication was performed with SamTools [32–34]. To reduce false-positive identification of ARGs, only gene accessions with reads aligning to > 80% of each reference nucleotide sequence were included in downstream analyses. Reads aligning to gene accessions requiring the confirmation of single nucleotide polymorphisms to confer resistance were excluded from downstream analysis.

Statistical analyses of resistome count data were performed using R v3.6.1 GUI 1.70 El Capitan build (7684) and RStudio v1.2.5001 [35, 36]. Custom scripting was used for processing of the data [37]. Major packages employed included metagenomeSeq [38, 39], tidyr [40], vegan [41], ggplot2 [42], tidyverse [43], phyloseq [44], data.table [45], and gridExtra [46]. MEGARes utilizes a hierarchical annotation structure with a label/level system to describe each gene accession by the “group” of genes it belongs to (based on sequence similarity, and not “gene” name), the “mechanism” by which it confers resistance, as well as the “class” of drug to which resistance is conferred. Therefore, references to specific groups in this study (e.g., TEM, SHV, ERMB, etc.) are denoting the label/level of “group”—representing sequence homology—and not gene names (e.g., tem, shv, ermB, etc.). Counts of alignments to resistance gene accessions were aggregated to the resistance drug “class,” “mechanism,” and “group” levels as per the MEGARes annotation structure to avoid irregularities with classifications at the “gene” level, particularly with the naming criteria for new resistance genes [47].

Richness and Inverse Simpson's diversity indices were calculated using vegan, and values were statistically compared using the Wilcoxon test function in R. Relative abundance (RA) 100% stacked barplots were created using mean counts at the resistance class and resistance mechanism levels for sample groups, with cut-offs of 1% RA. Cumulative sum scaling (CSS) normalized counts [38] were Hellinger-transformed, and Euclidean distances were calculated to create ordination plots with non-metric multidimensional scaling (NMDS). Analysis of Similarities (ANOSIM) was used to test for significant differences in ordination clustering of samples between timepoints [48]. The fitZig function (metagenomeSeq) was used to fit a zero-inflated Gaussian model and perform differential abundance testing [38, 39] of the various resistome features between timepoints. Sparsely present features (present in < 5% of samples) were removed [38] and the Benjamini–Hochberg correction was used to control False Discovery Rate [49]. Alpha = 0.05 was selected as the cut-off for statistical significance, and only q values were reported [50]. To account for false positive detection in low abundance resistome features, average expression values were determined (metagenomeSeq) and only features with an average expression value ≥ 1 were included at the class level and ≥ 3 at the mechanism level. Heatmaps were generated using CSS log₂ normalized counts, using the cut-offs for average expression and RA for at least one farm as inclusion criteria to visualize not only overall trends but also individual animal variation. Due to differences in farm management practices, changes in the farm management over time, and different AMD protocols used, statistical comparisons between Farms 1 and 2 in the Post period were not performed.

Several classes of ARGs were identified a priori as being of particular interest in relation to AMDs commonly used on dairy farms to treat diarrhea (macrolides, lincosamides, and streptogramins (MLS), sulfonamides, trimethoprim, aminoglycosides, and beta-lactams) or pneumonia in calves (MLS, beta-lactams, phenicols, and fluoroquinolones), or due to their critical importance in human medicine (beta-lactams, fluoroquinolones, and glycopeptides). Specific resistance mechanisms and groups belonging to these antimicrobial classes were also of a priori interest, including those relating to extended-spectrum beta-lactamases and other beta-lactamases (TEM, SHV, OXA, CTX, CMY, KPC, and NDM groups); the MCR group; or other similar colistin resistance; and all vancomycin and quinolone resistance mechanisms.

Results

Study Samples

A total of 31 fecal samples from dairy calves (1–30 days old) yielded 500 ng of DNA in ~100 μL per sample, which was the cut-off required for resistome target-enrichment and sequencing. Each sampling timepoint (Farm 1 Pre; Farm 1 Post; Farm 2 Pre; Farm 2 Post) had 7–8 samples.

Across all samples, sequencing generated 12 056 971 total raw reads, with a mean of 388 434 per sample (range: 56 160–557 284). Of these raw reads, a total of 10 489 045 paired-end reads passed screening. A mean of 338 356 (standard deviation = 105 338) and median of 372 370 (range: 51 850–491 070) reads classified as ARGs were present per fecal sample. Cumulatively, these reads aligned to 1192 ARG accessions, representing 15 AMR drug classes, 51 AMR mechanisms, and 212 AMR groups. Five AMD classes accounted for > 85% of sequences on both farms, across timepoints (multi-drug resistance (MDR), tetracyclines, beta-lactams, aminoglycosides, and sulfonamides; Figure 1). Similarly, five AMR mechanisms accounted for > 50% of sequences on both farms, and across timepoints (multi-drug efflux pumps, tetracycline resistance ribosomal protection proteins, class A beta-lactamases, sulfonamide-resistant dihydropteroate synthases, and multi-drug resistance regulator; Figure 2).

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Comparisons Between Farms (Pre)

Comparing Farms 1 and 2 at the Pre timepoint, there were distinct differences in RA across many resistance classes. Fourteen resistance classes were present on both farms, and 9/14 (64.3%) were statistically different between farms (Table 2). Of classes with an RA of > 1%, on Farm 1, there was an increased abundance of the MLS resistance class (q = 0.007). Conversely, on Farm 2, there were increased abundances of beta-lactam and fluoroquinolone resistance classes (q = 0.005 and q < 0.001, respectively).

TABLE 2 Comparison of Farm 1 and Farm 2 resistome mean antimicrobial resistance class relative abundance percentage (RA) and standard deviation (σ) at the Pre timepoint (n = 7–8) using a zero-inflated Gaussian model with differential abundance testing (* denoting statistical significance).

Class	Farm 1 Pre—RA % (σ)	Farm 2 Pre—RA % (σ)	q
MDR	18.2 (7.8)	24.1 (5.3)	0.04*
Tetracyclines	22.1 (6.4)	14.1 (6.5)	0.20
Aminoglycosides	20.0 (5.1)	16.4 (1.7)	0.81
Beta-lactams	14.4 (8.5)	26.7 (3.2)	0.005*
Sulfonamides	11.6 (1.8)	9.6 (2.6)	0.70
Phenicol	6.6 (4.3)	4.8 (1.6)	0.98
MLS	4.9 (2.5)	1.0 (2.8)	0.007*
CAP	1.0 (0.5)	1.1 (0.3)	0.13
Bacitracin	0.6 (0.4)	1.0 (0.3)	0.03*
Trimethoprim	0.5 (0.2)	0.2 (0.1)	0.02*
Fluoroquinolones	0 (0)	0.8 (0.0)	< 0.001*
Glycopeptides	0.1 (0.1)	0.0 (0.0)	0.004*

Comparisons Over Time Within Farm 1

The class-level ARG richness and diversity analyses for Farm 1 are depicted in Figures S1 and S2. There was a significant increase in class-level richness (p < 0.001) on Farm 1 in the Post period. No other timepoint comparisons or analyses were statistically significant (p > 0.05 for all other comparisons). Ordination using NMDS (Figure 3) to evaluate resistome composition at the AMR class level depicts clusters, separated by timepoint. The clusters for Farm 1 were distinct and separate from one another and were significant upon ANOSIM testing (R = 0.29, p = 0.02).

[IMAGE OMITTED. SEE PDF]

The RA of AMR classes for Farm 1 are shown in Figure 1. Pre and Post mean RA and differential abundance testing for all Farm 1 resistance classes are presented in Table 3. The tetracycline and fluoroquinolone resistance classes were differentially abundant and increased significantly between the Pre and Post timepoints (22.1%–27.4%, q = 0.02; 0%–0.1%, q < 0.001, respectively). No other classes were differentially abundant on Farm 1.

TABLE 3 Comparison of Farm 1 Pre- and Post-antimicrobial use intervention mean antimicrobial resistance class relative abundance percentage (RA) and standard deviation (σ; n = 7–8) using a zero-inflated Gaussian model with differential abundance testing (* denoting statistical significance).

Class	Pre—RA % (σ)	Post—RA % (σ)	q
MDR	18.2 (7.9)	20.4 (5.3)	0.08
Tetracyclines	22.1 (6.4)	27.4 (6.5)	0.02*
Aminoglycosides	20.0 (5.1)	19.2 (1.7)	0.08
Beta-lactams	14.4 (8.5)	6.4 (3.2)	0.32
Sulfonamides	11.6 (1.8)	11.1 (2.6)	0.13
Phenicol	6.6 (4.3)	6.7 (1.6)	0.10
MLS	4.9 (2.5)	6.6 (2.8)	0.08
CAP	1.0 (0.5)	0.8 (0.3)	0.36
Bacitracin	0.6 (0.4)	0.6 (0.3)	0.36
Trimethoprim	0.5 (0.2)	0.3 (0.1)	0.36
Fluoroquinolones	0 (0)	0.1 (0.0)	< 0.001*
Glycopeptides	0.1 (0.1)	0.1 (0.0)	0.18

Farm 1 resistance mechanisms with RAs ≥ 1% are shown in Figure 2. Seven resistance mechanisms were differentially abundant in the Post period, with increases seen for multi-drug efflux pumps (12.0%–14.9%, q = 0.03); tetracycline resistance ribosomal protection proteins (13.9%–17.6%, q = 0.02); tetracycline resistance MFS efflux pumps (7.4%–9.0%, q = 0.02); aminoglycoside O-phosphotransferases (6.7%–7.8%, q = 0.006); aminoglycoside N-acetyltransferases (1.1%–1.4%, q = 0.05); macrolide resistance efflux pumps (0.6%–2.0%, q = 0.006); and 23S rRNA methyltransferases (2.2%–3.6%, q = 0.05). Among resistance mechanisms meeting inclusion criteria, only five mechanisms numerically decreased in RA (class A and C beta-lactamases, multi-drug regulator, aminoglycoside O-nucleotidyltransferases, and macrolide phosphotransferases), and of these, only the macrolide phosphotransferases were statistically significant (1.3%–0.3%, q = 0.05). Pre and Post median RA and differential abundance testing for all Farm 1 resistance mechanisms are presented in Table S1.

Critical evaluation of the ARG mechanisms and groups belonging to the beta-lactam class of AMR was performed. On Farm 1, the RA of class D beta-lactamases increased from 0.0% to 0.3% (q = 0.04) in the Post period. The class D beta-lactamase group OXA was the only group present in the samples from Farm 1, and thus the increase at the mechanism level was solely driven by this group (3.7%–7.2% of all beta-lactams on Farm 1, q < 0.001). A numerical decrease in RA of class A beta-lactamases (10.4%–3.1%, q = 0.39) in the Post period was primarily driven by a statistically significant decrease in the RA of the TEM resistance group (47.1%–19.3% of all beta-lactams on Farm 1, q = 0.004). In contrast, the class A beta-lactamase group SHV had a numerically, but not statistically significant, increase in RA in the Post period (0%–0.4% of all beta-lactams on Farm 1, q = 0.06). While the class C beta-lactamases numerically decreased from 2.0% to 1.2% (q = 0.88) post-intervention, CMY was the only resistance group to have a statistically significant change on Farm 1 and decreased in RA in the Post period (7.3%–6.8% of all beta-lactams on Farm 1, q = 0.002). Pre and Post mean RA and differential abundance testing for all Farm 1 beta-lactam resistance groups are presented in Table 4.

TABLE 4 Comparison of Farm 1 Pre- and Post-antimicrobial use intervention mean antimicrobial resistance beta-lactamase group labels relative abundance percentage (RA), standard deviation (σ), median RA percentage, and range (%; n = 7–8) using a zero-inflated Gaussian model with differential abundance testing (* denoting statistical significance).

Mechanism	Group	Pre—Mean RA % (σ)	Pre—Median RA % (range)	Post—Mean RA % (σ)	Post—Median RA % (range)	q
Class A beta-lactamases	ACI	1.5 (3.2)	0.2 (0–8.7)	11.7 (8.6)	10.8 (0–27.4)	< 0.001*
Class A beta-lactamases	BLA1	0 (0)	0 (0–0)	3.5 (2.1)	3.7 (0.8–6.5)	< 0.001*
Class A beta-lactamases	BLAA	0 (0)	0 (0–0)	0.0 (0.1)	0 (0–0)	0.67
Class A beta-lactamases	BLAZ	0 (0)	0 (0–0)	0.4 (0.3)	0.3 (0.1–0.9)	0.27
Class A beta-lactamases	CARB	23.3 (40.7)	0 (0–96.2)	1.0 (0.6)	1.1 (0.2–1.7)	< 0.001*
Class A beta-lactamases	CBLA	0.0 (0.0)	0 (0–0.0)	0 (0)	0 (0–0)	0.24
Class A beta-lactamases	CEPA	1.4 (2.3)	0.0 (0–5.0)	0.1 (0.1)	0.1 (0.0–0.4)	0.08
Class A beta-lactamases	CFX	1.3 (3.1)	0 (0–8.2)	10.3 (11.6)	6.2 (0–25.8)	0.001*
Class A beta-lactamases	CTX	4.3 (3.8)	4.4 (0–10.4)	13.0 (10.5)	9.8 (2.1–29.4)	0.72
Class A beta-lactamases	ROB	0.1 (0.3)	0 (0–0.7)	4.7 (6.2)	0.8 (0–14.3)	0.004*
Class A beta-lactamases	SHV	0 (0)	0 (0–0)	0.4 (0.3)	0 (0–0.8)	0.06
Class A beta-lactamases	TEM	47.1 (37.8)	63.2 (0.1–88.0)	19.3 (20.3)	9.0 (0.1–51.3)	0.004 *
Class B beta-lactamases	L1	0 (0)	0 (0–0)	0.7 (0.5)	0.7 (0.1–1.4)	0.07
Class C beta-lactamases	ACT	0 (0)	0 (0–0)	0.0 (0.0)	0 (0–0.0)	0.90
Class C beta-lactamases	BLAEC	9.6 (18.3)	2.8 (0.7–51.0)	9.9 (4.9)	8.1 (4.9–18.5)	0.39
Class C beta-lactamases	CMY	7.3 (12.1)	0 (0–31.7)	6.8 (19.2)	0 (0–54.3)	0.002*
Class C beta-lactamases	SRT	0.4 (0.8)	0 (0–2.2)	0 (0)	0 (0–0)	0.10
Class D beta-lactamases	OXA	3.7 (9.0)	0.0 (0–24.1)	7.2 (5.2)	7.6 (0.2–15.1)	< 0.001*

The AMR class-level log₂ normalized count heatmap for Farm 1 is presented in Figure 4; the associated data are in Table S2. On Farm 1, most classes had similar median log₂ normalized count scores between timepoints. The range of log₂ scores was also narrower in the Post period.

[IMAGE OMITTED. SEE PDF]

The AMR mechanism-level log₂ normalized count heatmap for Farm 1 is presented in Figure 5. All AMR mechanism-level log₂ count data are presented in Table S3. On Farm 1, 48/51 total mechanisms were present; log₂ count scores were similar across most mechanisms at both timepoints, with a narrower range of scores in the Post period. Fourteen mechanisms had log₂ scores of zero in the Pre period, with positive scores in the Post period.

[IMAGE OMITTED. SEE PDF]

Comparisons Over Time Within Farm 2

The class-level ARG richness and diversity analyses for Farm 2 are depicted in Figures S1 and S2. No timepoint comparisons or analyses were found to be statistically significant (p > 0.05 for all comparisons). Ordination using NMDS to evaluate resistome composition at the AMR class level (Figure 3) depicts clusters, separated by timepoint. Despite the overlap of the Farm 2 clusters, as well as the larger distance distribution in the Post period, the Farm 2 clusters were statistically distinct upon ANOSIM testing (R = 0.20, p = 0.02).

The RA of AMR classes for Farm 2 are shown in Figure 1. Pre and Post mean RA and differential abundance testing for all Farm 2 resistance classes are presented in Table 5. The RA of the sulfonamide and fluoroquinolone resistance classes decreased in the Post period (9.6%–5.1%, q = 0.006; 0.8%–0.1%, q = 0.004, respectively). The beta-lactam resistance class also had a numerical decrease in RA in the Post period but failed to reach statistical significance (26.7%–18.1%, q = 0.06).

TABLE 5 Comparison of Farm 2 Pre- and Post-antimicrobial use intervention mean antimicrobial resistance class relative abundance percentage (RA) and standard deviation (σ; n = 8) using a zero-inflated Gaussian model with differential abundance testing (* denoting statistical significance).

Class	Pre—RA % (σ)	Post—RA % (σ)	q
MDR	24.1 (3.1)	23.9 (3.0)	0.58
Beta-lactams	26.7 (6.9)	18.1 (10.2)	0.06
Tetracyclines	14.1 (4.2)	23.0 (7.4)	0.42
Aminoglycosides	16.4 (1.6)	18.3 (4.8)	0.92
Sulfonamides	9.6 (2.5)	5.1 (3.6)	0.006*
Phenicol	4.8 (0.6)	5.8 (3.1)	0.56
CAP	1.1 (0.1)	1.0 (0.2)	0.49
MLS	1.0 (0.4)	2.6 (2.7)	0.92
Bacitracin	1.0 (0.2)	0.8 (0.2)	0.18
Fluoroquinolones	0.8 (0.3)	0.1 (0.1)	0.004*
Trimethoprim	0.2 (0.1)	0.3 (0.1)	0.92
Glycopeptides	0.0 (0.1)	0.0 (0.0)	0.11

Farm 2 resistance mechanisms RA ≥ 1% are shown in Figure 2. There were statistically significant increases in the Post period in the RA of the tetracycline resistance MFS efflux pumps (4.1%–6.9%, q = 0.03) and 23S rRNA methyltransferases (0.1%–0.8%, q = 0.04) resistance mechanisms. Decreases in the Post period in RA were noted in the class A beta-lactamases (21.6%–13.4%, q = 0.02), sulfonamide-resistant dihydropteroate synthases (9.9%–5.9%, q = 0.009), and quinolone resistance protein Qnr (0.8%–0.1%, q = 0.005) resistance mechanisms. Pre and Post median RA and differential abundance testing for all Farm 2 resistance mechanisms are presented in Table S4.

The decrease in RA in the Post period of class A beta-lactamases on Farm 2 was largely determined by decreases in the RA of the CTX and SHV resistance groups; the former approached significance, while the latter was statistically significant (32.6%–20.7% of all beta-lactams on Farm 2, q = 0.06; 0.6%–0% of all beta-lactams on Farm 2, q = 0.01, respectively). Pre and Post median RA and differential abundance testing for all Farm 2 beta-lactam resistance groups are presented in Table 6.

TABLE 6 Comparison of Farm 2 Pre- and Post-antimicrobial use intervention mean antimicrobial resistance beta-lactamase group labels relative abundance percentage (RA), standard deviation (σ), median RA percentage, and range (%; n = 8) using a zero-inflated Gaussian model with differential abundance testing (* denoting statistical significance).

Mechanism	Group	Pre—Mean RA % (σ)	Pre—Median RA % (range)	Post—Mean RA % (σ)	Post—Median RA % (range)	q
Class A beta-lactamases	ACI	0.1 (0.0)	0 (0–0.9)	0 (0)	0 (0–0)	0.02*
Class A beta-lactamases	BLA1	0 (0)	0 (0–0)	0.0 (0.0)	0 (0–0.0)	0.31
Class A beta-lactamases	CARB	0 (0)	0 (0–0)	1.7 (4.5)	0.2 (0–12.8)	0.05*
Class A beta-lactamases	CEPA	0.0 (0.0)	0.0 (0.0–0.1)	0.1 (0.2)	0.0 (0.0–0.6)	0.49
Class A beta-lactamases	CFX	1.0 (1.5)	0 (0–3.9)	3.1 (4.5)	0.7 (0–10.3)	0.93
Class A beta-lactamases	CTX	32.6 (17.0)	33.0 (14.4–59.1)	20.7 (13.3)	19.0 (2.2–44.3)	0.06
Class A beta-lactamases	ROB	0.0 (0.0)	0 (0–0.0)	4.4 (8.2)	0 (0–19.6)	< 0.001*
Class A beta-lactamases	SHV	0.6 (1.2)	0 (0–2.9)	0 (0)	0 (0–0)	0.01*
Class A beta-lactamases	TEM	49.8 (8.9)	50.3 (36.2–63.2)	38.9 (23.6)	45.3 (0.5–63.3)	0.13
Class B beta-lactamases	L1	0 (0)	0 (0–0)	0.2 (0.6)	0 (0–1.6)	0.05*
Class C beta-lactamases	ACT	0.0 (0.0)	0 (0–0.0)	0 (0)	0 (0–0)	0.02*
Class C beta-lactamases	BLAEC	2.5 (1.2)	2.1 (1.2–4.5)	5.7 (4.7)	2.1 (1.4–13.0)	0.5
Class C beta-lactamases	CMY	8.5 (10.4)	4.6 (0.1–25.6)	12.9 (14.7)	8.5 (0–33.1)	0.2
Class D beta-lactamases	OXA	0.0 (0.0)	0 (0–0.1)	0.0 (0.0)	0 (0–0.0)	0.26

The AMR class-level log₂ normalized count heatmap for Farm 2 is presented in Figure 4 and the associated data are in Table S5. Most classes in Farm 2 had similar median log₂ normalized count scores between timepoints. However, unlike Farm 1, Farm 2 had a narrower range of log₂ scores during the Pre period.

The AMR mechanism-level log₂ normalized count heatmap for Farm 2 is presented in Figure 5. All AMR mechanism-level log₂ count data are presented in Table S6. On Farm 2, 35/51 total mechanisms were present; log₂ count scores were similar across most mechanisms at both timepoints, with a narrower range of scores in the Pre period. Twenty-four out of the 35 mechanisms present (68.6%) showed resistance in 5/8 samples or greater across both timepoints.

Discussion

This study demonstrated there were disparate resistance patterns on each farm before AMU modification as well as within each farm, 1 year after the AMU protocol was implemented. The distinctive resistance patterns seen on each farm might reflect the diverse selection of AMDs used in the Pre period. Several studies have demonstrated that different AMDs have differing effects on the fecal microbiota and ARGs recovered [51, 52]. Alternatively, the varying management practices on each farm, or even different baseline microbiomes, could explain the disparities in the resistome seen between farms in the Pre period, as has been noted in several previous studies [7, 12, 53, 54]. The lack of concurrent controls, the limited number of farms, and the sampling of different animals between timepoints limit assessment of the effects of changes in AMU on AMR on these farms. Nevertheless, AMU—alongside other influencing factors—should be investigated to better understand the complexity and the potential role of antimicrobial stewardship approaches on AMR.

An association between age, AMR genes, and changes in microbial communities has been reported in calves, with lower detection rates of ARGs as calves age [8, 13, 14, 52, 55–57]. While all calves enrolled in the study were 1–30 days of age, precise age at sampling was not available. Significant changes in the fecal microbiota of calves between birth and 6 weeks of age (regardless of exposure to drug residues) have been previously reported [13], which encompasses the age of calves in our study. The microbiota also shifts both during the transition from feeding colostrum to feeding milk/milk replacer and during the first week of life [14]. Diet also affects resistome changes over the first few weeks of life, with differences seen in community structure during this timeframe [22]. If there were associations between the microbiota and resistome in the study subjects—independent of AMU—age and diet of the calves might have played a role in the types of ARGs recovered not only in the Pre period, but across collection times, for both farms.

The previous study of these farms found an increase in the bacterial richness in the Post period on both farms and that Farm 1 also had an increase in overall bacterial diversity in the Post period [15]. Several of the Post period differentially abundant taxa [15] have been associated with butyrate production [58, 59] and overall gut health in cattle [10]. Despite a purported shift towards a fecal microbiota associated with health, few differences in resistome RA were seen between timepoints on Farm 1. Numerical increases in RA of AMR were identified in the Post period for more than half the AMR classes. In fact, the total ARG counts recovered from the Farm 1 Post period were 1.2 times that of those recovered from the Pre period (~2.8 million counts vs. ~2.3 million counts). Statistically significant increases in the RA of fluoroquinolone and tetracycline resistance classes were observed in the Post period on Farm 1. These increases were also statistically significant for all fluoroquinolone mechanisms and for two-thirds of the tetracycline mechanisms. While an increase in the absolute counts of resistance recovered is incongruous with the large decrease in AMU on Farm 1, it could be related to the changes in bacterial richness and diversity and the microbial shift that were seen previously [15]. This microbial shift can potentially explain the change in the Farm 1 resistome composition as well as the increases in AMR seen on Farm 1 during the Post period. However, direct correlations cannot be made because the microbiome and resistome analyses were not run concurrently during DNA sequencing. Nonetheless, several other studies—across species—have noted the relationship between AMU and AMR is complex and that AMU or changes in AMU are not always reflected in the resulting ARGs and overall AMR [54, 60–63]. No distinction was made between intrinsic and acquired resistance in our study and so changes in the diversity of ARGs could be the result of the introduction of ARGs that are not responding to selection pressure whatsoever (e.g., efflux pumps). Co-selection and the continued use of some AMDs for treating diarrhea and/or pneumonia might also play a role in the continued persistence of some ARGs (e.g., ERMB and VANA) and could explain the disjointed AMU/AMR narrative [64]. Similarly, the use of AMDs on-farm for the treatment of cows (various conditions—AMD data not collected) could have contributed to the persistence of ARGs. However, it has also been shown that AMR persists in the environment with or without AMU [65, 66].

While changes in the resistome do occur over time on farms without any AMU intervention [22] as well as due to microbial shifts from AMU, Farm 2 appeared to have a reduction in ARGs in response to decreased AMU. Along with a 72% reduction in the number of calves treated for neonatal diarrhea due to the AMU intervention [15], a statistically significant reduction in the RA of the sulfonamide and fluoroquinolone resistance classes and a numerical decrease in both relative abundance and overall counts of beta-lactams was observed in the Post period. The change in protocol from using TMS and ceftiofur sodium in the Pre period to only using ceftiofur sodium as an alternative treatment some of the time in the Post period appears to be reflected in the data. Reductions for some mechanisms belonging to these classes were also observed and were statistically significant (class A beta-lactamases, sulfonamide-resistant dihydropteroate synthases, and quinolone resistance protein Qnr). These decreases are of interest because they demonstrate that specific reductions in AMU could potentially result in reductions of clinically important AMR classes, even in a complex environment where antimicrobials are used for other purposes in calves and cows. Similar reductions in AMR after reductions in AMU were also seen with the banning of avoparcin in Europe [67] and a reduction in 3rd-generation cephalosporin resistance in E. coli and Salmonella in chickens after restriction of prophylactic in ovo administration of ceftiofur [68]. While these studies have evaluated post-intervention assessments, there tends to be a narrow focus on a single type of bacteria and/or AMD, nor do they consider broader resistome changes. One defining strength of our work is that we evaluated changes to the full antibacterial resistome due to AMU modifications, which is uncommon in many studies evaluating AMR [20, 21, 69].

The concurrent and continued use of beta-lactam AMDs on-farm for calf pneumonia or in cows for other health conditions could have blunted the effect of reducing the use of ceftiofur to treat diarrhea. The presence of tetracycline ARGs is unsurprising given the extensive use of tetracycline AMDs in agriculture for over 70 years [70]. Many tetracycline resistance determinants are associated with plasmids—especially among gram-negative bacteria—which could partially explain their widespread distribution among bacterial species and explain some of the co-resistance seen with tetracyclines as well [71]. Tetracycline resistance is also common among several of the Gram-positive butyrate producers that were identified in the post period [15] and so the increase in this population of the microbiota might have brought along with it additional chromosomal ARGs. However, as has been reported previously, AMR has been found in various environments where human influence is not present [65, 66, 72, 73].

There are limitations to this study. One factor which hindered our ability to identify AMR changes related to the modifications in AMU treatment protocol was the relatively small sample size for each farm, at each timepoint, as well as the number of farms, and the lack of control farms (where no protocol change was implemented). Both individual animal and farm variation was noted in this study, and the smaller sample size could have affected the statistical power. It is possible that not all differences between timepoints were captured as significantly different in the statistics. Similarly, by only having two farms for comparison, broader trends in the resistome of dairy calves in southern Ontario might have been missed, particularly since inter-farm variation was clearly evident here and has been noted in previous work [23].

Additionally, different calves were sampled in the Pre and Post periods as the goal was to evaluate the farm-level effect of an intervention. Since sampling points were a year apart, changes in the resistome observed on both Farms 1 and 2 could have resulted from changes in bacterial populations rather than a true alteration of AMR. However, whether these changes were driven by resistance or microbiota shifts is to some degree irrelevant if they were associated with changes in AMU. There could have been differences in response to antimicrobial treatment based on changes in the farming practices and/or background microbiota. However, these were not quantified in this study, making it difficult to differentiate whether the changes in the resistome had more to do with changes in bacterial populations, or whether changes were due to individual animal differences, and thus a different set of AMR genes being present.

It is also possible that not enough time had elapsed from significant reductions in AMU to capture the resulting decreases in AMR. Given the protracted use of antimicrobials, not just on the farms in this study, but in agriculture at large, as well as several studies noting timeframes beyond 1 year to lose ARG carriership [21, 74, 75], more follow-up—at 2 and 5 years after the algorithm implementation—would have been useful to help determine the direction of AMR patterns and if further reductions could be seen with time.

Despite similar reductions in AMU on Farms 1 and 2, different resistance patterns were seen after implementing the decision-making algorithm. While our study design precluded being able to determine if the changes to the resistance patterns were a result of the decreased AMU or a shift in the microbiota, a result of the particular calves that were sampled, or due to changes in farm management practices, the trend on Farm 1 was towards increased abundance of several key AMR classes, while the trend on Farm 2 was towards decreased abundance of several key AMR classes. These diverging patterns of resistance highlight the complexity of AMR and that numerous factors can influence ARGs, making evaluation of changes to resistomes difficult. Although the presence of ARGs of a priori interest was noted, the abundance of these elements was low and likely attributable to a microbial shift and intrinsic resistance of the resulting bacteria, or from indirect effects of co-selection of ARGs. AMU policy and protocol modifications are becoming increasingly common, and studies such as this are crucial in monitoring the effects of changes to help maximize policy effectiveness.

Acknowledgments

The authors thank Cory Wolfe and Jake Castle from the Microbial Ecology Lab for performing the enriched library preparation and DNA sequencing.

Disclosure

Authors declare no off-label use of antimicrobials.

Ethics Statement

Approved by the University of Guelph Animal Care Committee: eAUP 3793. Authors declare human ethics approval was not needed.

Conflicts of Interest

The authors declare no conflicts of interest.

References

USDA, “Dairy 2007, Part I: Reference of Dairy Cattle Health and Management Practices in the United States United States Department of Agriculture Animal and Plant Health Inspection Service Veterinary Services (USDA‐APHIS‐VS),” 2005, 91–95.

United States Department of Agriculture, “Dairy 2014: Health and Management Practices on U.S. Dairy Operations,” 2014.

United States Department of Agriculture, “Dairy 2007—Heifer Calf Health and Management Practices on U.S. Operations,” 2010, USDA:APHIS:VS, CEAH Fort Collins, Co.

A. C. B. Berge, D. A. Moore, T. E. Besser, and W. M. Sischo, “Targeting Therapy to Minimize Antimicrobial Use in Preweaned Calves: Effects on Health, Growth, and Treatment Costs,” Journal of Dairy Science 92, no. 9 (2009): 4707–4714, https://doi.org/10.3168/jds.2009‐2199.

M. C. Windeyer, K. E. Leslie, S. M. Godden, D. C. Hodgins, K. D. Lissemore, and S. J. LeBlanc, “Factors Associated With Morbidity, Mortality, and Growth of Dairy Heifer Calves up to 3 Months of Age,” Preventive Veterinary Medicine 113 (2014): 231–240, https://doi.org/10.1016/j.prevetmed.2013.10.019.

D. E. Gomez, L. G. Arroyo, D. L. Renaud, L. Viel, and J. S. Weese, “A Multidisciplinary Approach to Reduce and Refine Antimicrobial Drugs Use for Diarrhoea in Dairy Calves,” Veterinary Journal 274 (2021): 105713, https://doi.org/10.1016/j.tvjl.2021.105713.

D. Gomez, L. Arroyo, M. Costa, L. Viel, and J. Weese, “Characterization of the Fecal Bacterial Microbiota of Healthy and Diarrheic Dairy Calves,” Journal of Veterinary Internal Medicine 31, no. 3 (2017): 928–939, https://doi.org/10.1111/jvim.14695.

Y. Uyeno, Y. Sekiguchi, and Y. Kamagata, “rRNA‐Based Analysis to Monitor Succession of Faecal Bacterial Communities in Holstein Calves,” Letters in Applied Microbiology 51, no. 5 (2010): 570–577, https://doi.org/10.1111/j.1472‐765X.2010.02937.x.

D. Klein‐Jöbstl, E. Schornsteiner, E. Mann, M. Wagner, M. Drillich, and S. Schmitz‐Esser, “Pyrosequencing Reveals Diverse Fecal Microbiota in Simmental Calves During Early Development,” Frontiers in Microbiology 5 (2014): 1–8, https://doi.org/10.3389/fmicb.2014.00622.

G. Oikonomou, A. G. V. Teixeira, C. Foditsch, M. L. Bicalho, V. S. Machado, and R. C. Bicalho, “Fecal Microbial Diversity in Pre‐Weaned Dairy Calves as Described by Pyrosequencing of Metagenomic 16S rDNA. Associations of Faecalibacterium Species With Health and Growth,” PLoS One 8, no. 4 (2013): e63157, https://doi.org/10.1371/journal.pone.0063157.

G. S. Slanzon, B. J. Ridenhour, D. A. Moore, et al., “Fecal Microbiome Profiles of Neonatal Dairy Calves With Varying Severities of Gastrointestinal Disease,” PLoS One 17 (2022): e0262317, https://doi.org/10.1371/journal.pone.0262317.

C. Foditsch, R. V. V. Pereira, J. D. Siler, C. Altier, and L. D. Warnick, “Effects of Treatment With Enrofloxacin or Tulathromycin on Fecal Microbiota Composition and Genetic Function of Dairy Calves,” PLoS One 14 (2019): e0219635, https://doi.org/10.1371/journal.pone.0219635.

R. V. V. Pereira, L. M. Carroll, S. Lima, et al., “Impacts of Feeding Preweaned Calves Milk Containing Drug Residues on the Functional Profile of the Fecal Microbiota,” Scientific Reports 8, no. 1 (2018): 554, https://doi.org/10.1038/s41598‐017‐19021‐2.

N. Amin and J. Seifert, “Dynamic Progression of the Calf's Microbiome and Its Influence on Host Health,” Computational and Structural Biotechnology Journal 19 (2021): 989–1001, https://doi.org/10.1016/j.csbj.2021.01.035.

D. E. Gomez, L. G. Arroyo, Z. Poljak, L. Viel, and J. S. Weese, “Implementation of an Algorithm for Selection of Antimicrobial Therapy for Diarrhoeic Calves: Impact on Antimicrobial Treatment Rates, Health and Faecal Microbiota,” Veterinary Journal 226 (2017): 15–25, https://doi.org/10.1016/j.tvjl.2017.06.009.

A. Garzon, R. Portillo, G. Habing, N. Silva‐del‐Rio, B. M. Karle, and R. V. Pereira, “Antimicrobial Resistance of Escherichia coli From Dairy Farms Participating in an Antimicrobial Stewardship Educational Program for Farm Employees,” Journal of Dairy Science 107, no. 3 (2024): 1645–1655, https://doi.org/10.3168/jds.2023‐23574.

R. M. Collis, P. J. Biggs, S. A. Burgess, A. C. Midwinter, G. Brightwell, and A. L. Cookson, “Impact of Systemic Antimicrobial Therapy on the Faecal Microbiome in Symptomatic Dairy Cows,” PLoS One 19 (2024): e0296290, https://doi.org/10.1371/journal.pone.0296290.

B. Awosile, J. McClure, J. Sanchez, J. C. Rodriguez‐Lecompte, G. Keefe, and L. C. Heider, “Salmonella Enterica and Extended‐Spectrum Cephalosporin‐Resistant Escherichia Coli Recovered From Holstein Dairy Calves From 8 Farms in New Brunswick, Canada,” Journal of Dairy Science 101, no. 4 (2018): 3271–3284, https://doi.org/10.3168/jds.2017‐13277.

A. Duse, K. P. Waller, U. Emanuelson, H. E. Unnerstad, Y. Persson, and B. Bengtsson, “Risk Factors for Antimicrobial Resistance in Fecal Escherichia coli From Preweaned Dairy Calves,” Journal of Dairy Science 98, no. 1 (2015): 500–516, https://doi.org/10.3168/jds.2014‐8432.

J. A. Afema, M. A. Davis, and W. M. Sischo, “Antimicrobial Use Policy Change in Pre‐Weaned Dairy Calves and Its Impact on Antimicrobial Resistance in Commensal Escherichia Coli: A Cross Sectional and Ecological Study,” BMC Microbiology 19 (2019): 217, https://doi.org/10.1186/s12866‐019‐1576‐6.

M. de Lagarde, J. M. Fairbrother, M. Archambault, et al., “Impact of a Regulation Restricting Critical Antimicrobial Usage on Prevalence of Antimicrobial Resistance in Escherichia coli Isolates From Fecal and Manure Pit Samples on Dairy Farms in Québec, Canada,” Frontiers in Veterinary Science 9 (2022): 9, https://doi.org/10.3389/fvets.2022.838498.

J. Liu, D. H. Taft, M. X. Maldonado‐Gomez, et al., “The Fecal Resistome of Dairy Cattle Is Associated With Diet During Nursing,” Nature Communications 10, no. 1 (2019): 4406, https://doi.org/10.1038/s41467‐019‐12111‐x.

B. J. Haley, S. W. Kim, S. Salaheen, E. Hovingh, and J. A. S. Van Kessel, “Differences in the Microbial Community and Resistome Structures of Feces From Preweaned Calves and Lactating Dairy Cows in Commercial Dairy Herds,” Foodborne Pathogens and Disease 17, no. 8 (2020): 494–503, https://doi.org/10.1089/fpd.2019.2768.

M. J. Drake, S. G. Daniel, L. D. Baker, et al., “Effects of Dietary Zinc on the Gut Microbiome and Resistome of the Gestating Cow and Neonatal Calf,” Animal Microbiome 6, no. 1 (2024): 39, https://doi.org/10.1186/s42523‐024‐00326‐3.

N. R. Noyes, M. E. Weinroth, J. K. Parker, et al., “Enrichment Allows Identification of Diverse, Rare Elements in Metagenomic Resistome‐Virulome Sequencing,” Microbiome 5, no. 1 (2017): 142, https://doi.org/10.1186/s40168‐017‐0361‐8.

Canadian Veterinary Medical Association, “Access the CVMA Veterinary Guidelines for AMU on the Firstline App,” 2022, Firstline.

Canadian Veterinary Medical Association, “Diarrheic Beef Calf (<30 Days Old) – Antimicrobial Treatment Decision Tree,” 2022, CVMA, https://www.canadianveterinarians.net/veterinary‐resources/pugs/resource‐library/#3‐prescription.

S. M. Lakin, C. Dean, N. R. Noyes, et al., “MEGARes: An Antimicrobial Resistance Database for High Throughput Sequencing,” Nucleic Acids Research 45, no. D1 (2017): D574–D580, https://doi.org/10.1093/nar/gkw1009.

C. J. Dean, “Welcome to AmrPlusPlus's Documentation!,” 2016, http://megares.meglab.org/amrplusplus/latest/html/index.html.

A. M. Bolger, M. Lohse, and B. Usadel, “Trimmomatic: A Flexible Trimmer for Illumina Sequence Data,” Bioinformatics 30, no. 15 (2014): 2114–2120, https://doi.org/10.1093/bioinformatics/btu170.

S. Merchant, D. E. Wood, and S. L. Salzberg, “Unexpected Cross‐Species Contamination in Genome Sequencing Projects,” PeerJ 2 (2014): e675, https://doi.org/10.7717/peerj.675.

H. Li, “Aligning Sequence Reads, Clone Sequences and Assembly Contigs With BWA‐MEM,” arXiv:13033997 [q‐bioGN], 2013, 1–3, http://arxiv.org/abs/1303.3997.

H. Li, B. Handsaker, A. Wysoker, et al., “The Sequence Alignment/Map Format and SAMtools,” Bioinformatics 25, no. 16 (2009): 2078–2079, https://doi.org/10.1093/bioinformatics/btp352.

H. Li and R. Durbin, “Fast and Accurate Long‐Read Alignment With Burrows‐Wheeler Transform,” Bioinformatics 26, no. 5 (2010): 589–595, https://doi.org/10.1093/bioinformatics/btp698.

R Core Team, “R: A Language and Environment for Statistical Computing,” 2019, https://www.R‐project.org/.

RStudio Team, “RStudio: Integrated Development for R,” 2019, http://www.rstudio.com/.

E. Doster, “Microbial‐Ecology‐Group/Guelph_Calf_Diarrhea_Metagenomic_Analyses,” Microbial Ecology Group GitHub, 2020, https://github.com/Microbial‐Ecology‐Group/Guelph_calf_diarrhea_metagenomic_analyses.

J. N. Paulson, C. O. Stine, H. C. Bravo, and M. Pop, “Differential Abundance Analysis for Microbial Marker‐Gene Surveys,” Nature Methods 10, no. 2 (2013): 1200–1202, https://doi.org/10.1038/nmeth.2658.

J. N. Paulson, H. Talukder, M. Pop, and H. C. Bravo, “metagenomeSeq: Statistical Analysis for Sparse High‐Throughput Sequencing. Bioconductor Package: 1280,” 2013, http://cbcb.umd.edu/software/metagenomeSeq.

H. Wickham and L. Henry, “tidyr: Tidy Messy Data,” R Package Version 100, 2019, https://cran.r‐project.org/package=tidyr.

J. Oksanen, G. L. Simpson, P. Solymos, et al., “vegan: Community Ecology Package,” R Package Version 256, 2019, https://cran.r‐project.org/package=vegan.

H. Wickham, “ggplot2: Elegant Graphics for Data Analysis,” 2016, https://ggplot2.tidyverse.org.

H. Wickham, “tidyverse: Easily Install and Load the ‘Tidyverse’,” R Package Version 121, 2017, https://CRAN.R‐project.org/package=tidyverse.

P. J. McMurdie and S. Holmes, “Phyloseq: An R Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Data,” PLoS One 8, no. 4 (2013): e61217, https://doi.org/10.1371/journal.pone.0061217.

M. Dowle and A. Srinivasan, “data.table: Extension of ‘data.frame’,” R Package Version 1126, 2019, https://CRAN.R‐project.org/package=data.table.

B. Auguie, “gridExtra: Miscellaneous Functions for ‘Grid’ Graphics,” R Package Version 23, 2017, https://CRAN.R‐project.org/package=gridExtra.

R. M. Hall and S. Schwarz, “Resistance Gene Naming and Numbering: Is It a New Gene or Not?,” Journal of Antimicrobial Chemotherapy 71, no. 3 (2016): 569–571, https://doi.org/10.1093/jac/dkv351.

K. R. Clarke, “Non‐Parametric Multivariate Analyses of Changes in Community Structure,” Australian Journal of Ecology 18, no. 1 (1993): 117–143, https://doi.org/10.1111/j.1442‐9993.1993.tb00438.x.

Y. Benjamini and Y. Hochberg, “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing,” Journal of the Royal Statistical Society 57, no. 1 (1995): 289–300, https://doi.org/10.2307/2346101.

J. D. Storey, “A Direct Approach to False Discovery Rates,” Journal of the Royal Statistical Society, Series B: Statistical Methodology 64, no. 3 (2002): 479–498, https://doi.org/10.1111/1467‐9868.00346.

A. M. R. Grønvold, Y. Mao, T. M. L'Abée‐Lund, et al., “Fecal Microbiota of Calves in the Clinical Setting: Effect of Penicillin Treatment,” Veterinary Microbiology 153, no. 3–4 (2011): 354–360, https://doi.org/10.1016/j.vetmic.2011.05.040.

R. V. Pereira, J. D. Siler, J. C. Ng, M. A. Davis, Y. T. Grohn, and L. D. Warnick, “Effect of On‐Farm Use of Antimicrobial Drugs on Resistance in Fecal Escherichia coli of Preweaned Dairy Calves,” Journal of Dairy Science 97, no. 12 (2014): 7644–7654, https://doi.org/10.3168/jds.2014‐8521.

J. S. Weese and M. Jelinski, “Assessment of the Fecal Microbiota in Beef Calves,” Journal of Veterinary Internal Medicine 31, no. 1 (2017): 176–185, https://doi.org/10.1111/jvim.14611.

E. Doster, L. J. Pinnell, N. R. Noyes, et al., “Evaluating the Effects of Antimicrobial Drug Use on the Ecology of Antimicrobial Resistance and Microbial Community Structure in Beef Feedlot Cattle,” Frontiers in Microbiology 13 (2022): 13, https://doi.org/10.3389/fmicb.2022.970358.

M. Rey, F. Enjalbert, S. Combes, L. Cauquil, O. Bouchez, and V. Monteils, “Establishment of Ruminal Bacterial Community in Dairy Calves From Birth to Weaning Is Sequential,” Journal of Applied Microbiology 116, no. 2 (2014): 245–257, https://doi.org/10.1111/jam.12405.

M. L. Hennessy, N. Indugu, B. Vecchiarelli, et al., “Temporal Changes in the Fecal Bacterial Community in Holstein Dairy Calves From Birth Through the Transition to a Solid Diet,” PLoS One 15 (2020): 1–16, https://doi.org/10.1371/journal.pone.0238882.

N. R. Noyes, X. Yang, L. M. Linke, et al., “Characterization of the Resistome in Manure, Soil and Wastewater From Dairy and Beef Production Systems,” Scientific Reports 6 (2016): 6, https://doi.org/10.1038/srep24645.

P. Górka, Z. M. Kowalski, R. Zabielski, and P. Guilloteau, “Invited Review: Use of Butyrate to Promote Gastrointestinal Tract Development in Calves,” Journal of Dairy Science 101, no. 6 (2018): 4785–4800, https://doi.org/10.3168/jds.2017‐14086.

A. Bedford and J. Gong, “Implications of Butyrate and Its Derivatives for Gut Health and Animal Production,” Animal Nutrition 4, no. 2 (2018): 151–159, https://doi.org/10.1016/j.aninu.2017.08.010.

P. Rovira, T. McAllister, S. M. Lakin, et al., “Characterization of the Microbial Resistome in Conventional and ‘Raised Without Antibiotics’ Beef and Dairy Production Systems,” Frontiers in Microbiology 10 (2019): 1–12, https://doi.org/10.3389/fmicb.2019.01980.

G. Van Den Bunt, A. C. Fluit, M. P. Spaninks, et al., “Faecal Carriage, Risk Factors, Acquisition and Persistence of ESBL‐Producing Enterobacteriaceae in Dogs and Cats and Co‐Carriage With Humans Belonging to the Same Household,” Journal of Antimicrobial Chemotherapy 75, no. 2 (2020): 342–350, https://doi.org/10.1093/jac/dkz462.

V. Dazio, A. Nigg, J. S. Schmidt, et al., “Duration of Carriage of Multidrug‐Resistant Bacteria in Dogs and Cats in Veterinary Care and Co‐Carriage With Their Owners,” One Health 13 (2021): 100322, https://doi.org/10.1016/j.onehlt.2021.100322.

S. Rahman and A. Hollis, “The Effect of Antibiotic Usage on Resistance in Humans and Food‐Producing Animals: A Longitudinal, One Health Analysis Using European Data,” Frontiers in Public Health 11 (2023): 1170426, https://doi.org/10.3389/fpubh.2023.1170426.

F. M. Aarestrup, “Characterization of Glycopeptide‐Resistant Enterococcus Faecium (GRE) From Broilers and Pigs in Denmark: Genetic Evidence That Persistence of GRE in Pig Herds Is Associated With Coselection by Resistance to Macrolides,” Journal of Clinical Microbiology 38 (2000): 2774–2777.

V. M. D'Costa, C. E. King, L. Kalan, et al., “Antibiotic Resistance Is Ancient,” Nature 477, no. 7365 (2011): 457–461, https://doi.org/10.1038/nature10388.

L. J. Pinnell, G. Kupier, K. L. Huebner, et al., “More Than an Anthropogenic Phenomenon: Antimicrobial Resistance in Ungulates From Natural and Agricultural Environments,” Science of the Total Environment 858 (2023): 159789, https://doi.org/10.1016/j.scitotenv.2022.159789.

DANMAP, DANMAP 2013 ‐ Use of Antimicrobial Agents and Occurrence of Antimicrobial Resistance in Bacteria From Food Animals, Food and Humans in Denmark (Statens Serum Institut, 2013).

Public Health Agency of Canada, “Salmonella Heidelberg – Ceftiofur‐Related Resistance in Human and Retail Chicken Isolates,” 2006.

K. L. Tang, N. P. Caffrey, D. B. Nóbrega, et al., “Comparison of Different Approaches to Antibiotic Restriction in Food‐Producing Animals: Stratified Results From a Systematic Review and Meta‐Analysis,” BMJ Global Health 4, no. 4 (2019): e001710, https://doi.org/10.1136/bmjgh‐2019‐001710.

M. M. Swann, K. L. Blaxter, and H. I. Field, “Report to the Joint Committee on the Use of Antibiotics in Animal Husbandry and Veterinary Medicine,” 1969.

M. C. Roberts, “Tetracycline Resistance Determinants: Mechanisms of Action, Regulation of Expression, Genetic Mobility, and Distribution,” FEMS Microbiology Reviews 19 (1996): 1–24, https://doi.org/10.1016/0168‐6445(96)00021‐6.

K. Bhullar, N. Waglechner, A. Pawlowski, et al., “Antibiotic Resistance Is Prevalent in an Isolated Cave Microbiome,” PLoS One 7, no. 4 (2012): 1–11, https://doi.org/10.1371/journal.pone.0034953.

J. Perry, N. Waglechner, and G. Wright, “The Prehistory of Antibiotic Resistance,” Cold Spring Harbor Perspectives in Medicine 6, no. 6 (2016): 1–8, https://doi.org/10.1101/cshperspect.a025197.

P. F. M. Teunis, E. G. Evers, P. D. Hengeveld, C. M. Dierikx, C. C. C. H. Wielders, and E. Van Duijkeren, “Time to Acquire and Lose Carriership of ESBL/ pAmpC Producing E. coli in Humans in The Netherlands,” PLoS One 13, no. 3 (2018): e0193834, https://doi.org/10.1371/journal.pone.0193834.

R. A. Horton, D. Duncan, L. P. Randall, et al., “Longitudinal Study of CTX‐M ESBL‐Producing E. coli Strains on a UK Dairy Farm,” Research in Veterinary Science 109 (2016): 107–113, https://doi.org/10.1016/j.rvsc.2016.09.018.