We tested frozen swab heads from deep nasopharyngeal (DNP) swabs collected from 20 beef calves stored at −80°C as part of a larger sample collection. The samples used in this experiment were selected for a range of culture and phenotypic AST outcomes. Overall, 10 of the 20 samples were collected in 2020 from a single pen of calves sampled at 6 days on feed (DOF). Six samples were collected from calves sampled at 13 DOF from other pens from the same study, where calves arrived from October to early December 2020. Four samples were collected in 2021 from different early-, mid-, and late-filled feedlot pens (Supplementary Table S1A). The research protocols and procedures for this study were approved by the University of Saskatchewan Animal Care Committee (AUP 20190069).

Samples were collected at a research feedlot operated by the University of Saskatchewan from 1,600 recently weaned mixed-breed steers purchased from a western Canadian auction market in the fall of 2020 and 2021. Calves were restrained in a hydraulic chute with a neck-extender, and three DNP swabs were collected from each calf from alternating nostrils (collection protocol details in Supplementary material). All three DNP swabs per calf were pooled in the same vial.

Samples were transported to the University of Saskatchewan for processing. The pooled samples (three swabs from each calf) were vortexed for 1 min to release biomass from the swab to the transport medium, and an aliquot was submitted to Prairie Diagnostic Services, Inc. (PDS; Saskatoon, Saskatchewan, Canada) for culture and AST.

In total, 10 μL inoculation loop of sample was cultured on 5% Columbia sheep blood (BA) and chocolate agar (CHOC) plates and incubated at 35°C for 18 h in an environment of 5% CO₂ to isolate M. haemolytica, P. multocida, and H. somni. Bacterial colonies were incubated for 18 h and 42 h. Isolates of interest were identified using the MALDI-TOF MS Microflex LT instrument (Bruker Daltonik, Bremen, Germany) and the MALDI-TOF MS Biotyper Microflex LT Compass version 1.4 software with MSP library, according to the manufacturer’s guidelines. Isolate abundance was scored using a semi-quantitative scale (trace, 1+ to 4+) based on counts of visible colonies following streaking plates by quadrants (identification and quality control details in Supplementary material).

Susceptibility was measured by a commercially available serial broth microdilution panel using the Sensititre^™ platform and the commercially available bovine BOPO7F Plate (ThermoFisher Scientific^™, Thermo Fisher Scientific, Waltham, Massachusetts, United States), following the manufacturer’s instructions and recommendations for testing and quality control (antimicrobials and concentration ranges in Supplementary material). E. coli ATCC 25922, Staphylococcus aureus ATCC 29213, and H. somni ATCC 700025 were used to evaluate the performance. The minimum inhibitory concentration (MIC) plate was placed and read on the BIOMIC^® V3 microplate reader. The MIC value was considered equal to the lowest concentration of antimicrobial that inhibited visible growth.

The MICs were compared against breakpoints designated by the Clinical and Laboratory Standards Institute (CLSI) for the pathogens of interest. Isolates with MIC values considered intermediate were categorized as “susceptible” for all analyses (CLSI, 2023).

For samples collected in 2020, 2 mL of the remaining transport medium was centrifuged at 4,000 × g for 10 min to pellet biomass. Then, 900 μL of supernatant was decanted, and the pellet was resuspended in the remaining transport medium and host-depleted using the HostZERO™ Microbial DNA Kit (Zymo Research, Irvine, California, United States), according to the manufacturer’s recommended protocol and extracted using the QIAGEN DNEasy Blood and Tissue Vacuum Kit (QIAGEN, Hilden, Germany). The dry, swab heads were then stored in cryovials at −80°C without media or cryoprotectant.

Dry swab heads for all samples/years were removed from −80°C for enrichment. The swab head was thawed briefly at room temperature and placed aseptically into sterile glass vials with 14 mL of BHI medium containing 1% glucose and a stir bar. Glass vials were sealed with air permeable Air ‘o Top membranes (Thomson Instrument Company, Oceanside, California, United States). Cultures were grown with vigorous aeration at 35°C for 22 h. During this period, 1 mL of each sample was retrieved at 0 h and 8 h and, subsequently, every 2 h until the sample had reached the stationary phase of growth. Stationary phase was determined by optical density (OD) that indicated culture saturation at three consecutive timepoints. Each 1 mL sample was pelleted at 4,000 × g for 10 min, and the medium was discarded. Bacterial pellets were stored at −20°C overnight. Nucleic acid extraction of the resulting enriched bacterial pellets was performed using the Gentra Puregene Buccal Cell Kit (QIAGEN, Hilden, Germany), according to the manufacturer’s instructions.

DNA concentration was determined using Qubit™ 1x dsDNA Broad Range (BR) or High Sensitivity (HS) Assay Kits (Invitrogen, Carlsbad, California, United States), according to the manufacturer’s specifications. The extracted DNA was stored at 4°C for <1 month in TE-buffer until library preparation and qPCR.

To identify M. haemolytica, P. multocida, and H. somni, qPCR was performed in triplicate on the Aria MX (Agilent Technologies, Santa Clara, California, United States) using Taqman^™ Fast Advanced MasterMix (Invitrogen, Carlsbad, California, United States). All samples were normalized to 10 ng/μl for amplification and compared against a standard curve. Total microbial load was quantified by targeting the bacterial 16S rRNA gene (Nadkarni et al., 2002), and the abundance of organisms of interest was quantified by targeting species-specific marker genes (Kishimoto et al., 2017). The results from qPCR were examined to identify time point(s) with the highest concentrations among all BRD organisms of interest for comparative testing after the stationary phase of growth was attained.

Two candidate optimal enrichment timepoints were identified based on qPCR results for M. haemolytica, P. multocida, and H. somni. DNA from the optimal timepoints and an aliquot of DNA that was extracted without enrichment were selected for further testing. Sequencing libraries were prepared using the ONT ligation kit SQK-LSK109 and native barcoding kit (EXP-NBD104 and EXP-NBD114), according to the manufacturer’s instructions, with the following reductions to reaction volumes: repair and end prep reactions were scaled to 15 μL, and barcode ligation was scaled to 20 μL. Eight barcoded samples were normalized and pooled into each library, resulting in seven sequencing runs, each run containing eight samples. Sequencing libraries were quantified with the Qubit^™ 1x High Sensitivity (HS) Assay Kit (Invitrogen, Carlsbad, California, United States). Overall, 200 ng of each prepared library was loaded onto an FLO-MIN106 flow cell and sequenced on an ONT GridION device for 72 h.

ONT GridION default run parameters and high-accuracy Guppy (v4.0) real-time base calling were used to process raw signal data and remove reads with an average quality < Q7, after which terminal and internal adapters in split reads were removed with Porechop v0.2.4 (Wick et al., 2017). Reads shorter than 100 bp were removed using NanoFilt v2.6.0, and sequence statistics were calculated using NanoStat v1.5.0 (De Coster et al., 2018).

Kraken2 v2.0.8-beta (Wood et al., 2019) was used to classify host and non-host reads with a custom database. The database included all complete genomes in NCBI RefSeq for the bacterial, viral, and archaeal domains on 17 October 2020, as well as the Bos taurus reference genome assembly ARS-UCD1.2_Btau5.0.1Y, which consists of the ARS-UCD1.2 genome assembly (Elsik et al., 2009; Rosen et al., 2020). Following classification, reads were divided into two groups using the KrakenTools v1.0 utility extract_kraken_reads.py: those assigned to the B. taurus genome were placed into the “host” dataset and those were classified as any other taxa or were unclassified and placed into the “non-host” dataset. Chimeric reads were retrieved from the host dataset as described by Freeman et al. (2022). Taxonomic abundance estimates for organisms of interest were computed by Bracken v2.5 (Lu et al., 2017) from the “non-host” dataset after removing host-like sequences using the KrakenTools script filter_bracken_out.py. The number of reads and total base pairs were reported for BRD bacteria, such as M. haemolytica, P. multocida, H. somni, M. bovis, and B. trehalosi.

Antimicrobial resistance genes were identified in non-host reads using ABRicate v.1.0.1 (Seemann, n.d.) and AMRFinderPlus v3.9.8 (Feldgarden et al., 2019), both with the NCBI Bacterial Antimicrobial Resistance Reference Gene Database (PRJNA313047, version 2020-12-17). ABRicate was also run using the Comprehensive Antimicrobial Resistance Database (CARD) (Alcock et al., 2019) and MEGARes 2.0 database (Doster et al., 2019). For AMRFinderPlus, the minimum percent identity and percent coverage thresholds were set to 60%, as it is more stringent in reporting ARGs, and the -plus option was used to direct the program to search for genes involved in virulence, biocide, heat, metal, and acid resistance. Default parameters were used for ABRicate (80% minimum percent identity and percent coverage).

Theoretical genome coverage was calculated as the sum of the lengths of reads classified as a particular BRD pathogen divided by the size of its reference genome (M. haemolytica: 2.8 Mb [NCBI GCF_002285575.1], P. multocida: 2.3 Mb [NCBI GCF_002073255.2], H. somni: 2.2 Mb [NCBI GCF_000019405.1], M. bovis: 0.9 Mb [NCBI GCF_001930225.1], and B. trehalosi: 2.3 Mb [NCBI GCF_000521725.1]).

Sequence statistics (number of reads, total base pairs, and theoretical genome coverage) and qPCR results for organism detection were summarized as medians. qPCR copy numbers and sequence statistics were compared between culture-positive and culture-negative samples for no enrichment and 10h and 14h enrichment using the Wilcoxon rank sum test for each organism (StataSE ver 18.0, StataCorp, College Station, TX). The Wilcoxon signed-rank test was used to compare the total base pairs detected for each organism among matched pairs of samples that were unenriched to samples enriched for 10h and 14h. The number of total base pairs was also compared with qPCR concentration for the 10h and 14h enrichment scenarios using Spearman’s correlation coefficient.

Concordance was assessed between genomic detection, defined by surpassing a threshold of sequencing reads per organism, and the traditional culture results (either positive or negative) using the kappa statistic for each organism and enrichment duration. In the present study, the cutoff representing species detection with metagenomics was estimated based on the distribution of read counts to optimize the distinction between culture positive and negative samples, leveraging previously reported metagenomic sequencing examples from the literature that used cutoffs of 100 and 1,000 reads (Zhang et al., 2022; Liu et al., 2023). Kappa was interpreted as 0.81–1.0 almost perfect agreement; 0.61–0.80 substantial agreement, 0.40–0.60 moderate agreement, 0.21–0.40 fair agreement, and 0.01–0.20 none to slight agreement (Dohoo et al., 2009).

For AMR and ARG detection for antimicrobials of interest in managing BRD, statistical analyses included kappa to assess agreement between MIC-based AST results and detection of ARGs, Wilcoxon rank sum test to assess differences in ARG numbers between isolates with and without phenotypic AMR, and the Wilcoxon singed-rank test to assess differences in ARGs among unenriched samples and samples enriched for 10h and 14 h. p-values ≤ 0.05 were considered significant.

Sequential qPCR testing was used initially to screen samples from all time periods and select candidate enrichment time points for metagenomic sequencing. Total bacterial abundance continued to increase up to 10 h of incubation, as determined by qPCR targeting the 16S rRNA gene (Figures 1A,B). M. haemolytica plateaued after 8 h. P. multocida and H. somni were highest at 14 h. Samples from 10 and 14 h were selected for metagenomic sequencing; 10 h was the first time point after total bacterial abundance peaked and 14 h represented as the highest joint concentrations of M. haemolytica, P. multocida, and H. somni.

After removal of short (<100 bp) and low-quality (<Q7) sequences, the median total base pairs of data available for analyses, number of total reads, read lengths, and read quality were summarized for unenriched, 10- and 14-h enriched samples (Supplementary Tables S1A,B). The median read lengths were longer for enriched than unenriched samples. Non-selective enrichment increased the percentage of non-host DNA in most samples (Table 1). The percentage of non-host sequence measured by total base pairs and total non-host reads was less than 7% for all non-enriched samples. The median percentage of non-host sequence was higher for enriched samples; 54% (range 7–93%) for 10 h and 61% (range 38–94%) for 14 h (Table 1).

Bacterial enrichment increased the percentage of sequence data from BRD pathogens relative to unenriched samples. The increase was two to almost four orders of magnitude higher for read numbers and total sequence length (bp) for M. haemolytica and two orders of magnitude for P. multocida (Table 2, Figure 2, Supplementary Tables S2A,B). There were no comparable culture-positive unenriched samples for H. somni (Supplementary Table S2C).

In the absence of enrichment, only one sample produced >1x theoretical genome coverage for any organism (Supplementary Tables S2A–C). This sample had 2.8x coverage for M. haemolytica (Supplementary Table S2A) and was the only sample where a culture was semi-quantitatively scored as 3+ out of a possible 4+ abundance of colony growth on the culture plate (Table 2). In contrast, at 14 h, M. haemolytica was detected at >30x theoretical coverage in 11 of 19 culture-positive samples (all samples >1.8x coverage) (Supplementary Table S2A). At 14 h, P. multocida was detected at >1.6x theoretical coverage in five of seven culture-positive samples (Supplementary Table S2B), and H. somni was detected at >8x theoretical coverage in all four culture-positive samples (Supplementary Table S2C). The median theoretical coverage was ≤0.1 for culture negative samples of P. multocida and H. somni at 10h and 14 h (Supplementary Tables S2B,C). There was only one sample culture negative for M. haemolytica with a theoretical coverage of 1.3 for 10 h and 1.9 for 14 h (Supplementary Table S2A). Notably, this sample was also qPCR positive.

M. bovis was only detected in enriched samples (6/20 at 10 h and 2/20 at 14 h) and B. trehalosi was detected in more enriched (12/20 at 10 h and 18/20 at 14 h) than unenriched samples (3/16) (p < 0.05) (additional details are included in Supplementary material).

M. haemolytica had been cultured from 19 of 20 samples, followed by P. multocida (n = 7) and H. somni (n = 4). Two BRD pathogens were co-isolated from 10 of 20 samples (Supplementary Tables S2A–C).

Culture detection for each of the three organisms of interest was compared with taxonomic classification of sequencing data (Table 2). All 16 unenriched samples (Table 2, Supplementary Table S2A) were M. haemolytica positive, and kappa could not be calculated. One unenriched sample with a 3+ abundance culture score had >1,000 M. haemolytica reads and the rest had <250 reads. In total, 6 of the 16 culture-positive unenriched samples from 2020 had <10 reads for M. haemolytica.

As there was only one culture-negative sample, the estimation of kappa for M. haemolytica was limited for the 20 enriched samples. However, there was substantial agreement (κ = 0.64) between culture positivity and samples with >1,000 reads for 14h. There was very good correlation between the qPCR results and total base pairs sequenced for M. haemolytica at 10 h (ρ = 0.86) and 14h (ρ = 0.87) (Table 2, Supplementary Table S2A).

There were no samples with >100 P. multocida reads in the unenriched group (Table 2, Supplementary Table S2B). If samples were considered positive at ≥1 P. multocida reads, agreement of unenriched sequencing with culture results was poor (κ = 0.11). Agreement between detection of >100 reads and culture positivity was substantial (κ = 0.68) at 14 h and moderate (κ = 0.56) at 10 h (Table 2). At 14 h, two P. multocida culture-negative samples had >100 reads; one was strong positive for P. multocida on qPCR suggesting a false-negative culture (Supplementary Table S2B). At 10h and 14h, culture-positive-enriched samples had higher (p < 0.01) numbers of P. multocida reads, total base pairs, and theoretical genome coverage than culture-negative samples; whereas, there were no significant differences without enrichment (Table 2). There was also good correlation between qPCR results and total base pairs for P. multocida at 10h (ρ = 0.77) and 14 h (ρ = 0.79) (Table 2).

For H. somni (Table 2, Supplementary Table S2C), all 2020 samples were culture-negative. In the 10h- and 14h-enriched samples, the agreement between total H. somni reads >100 and culture-positive samples was substantial (κ = 0.63, κ = 0.74). Correlation was satisfactory to good between qPCR results and total base pairs of H. somni at 10 (ρ = 0.49) and 14 h (ρ = 0.68) (Table 2). The 10h- and 14h-enriched samples had higher (p < 0.001) numbers of H. somni reads, total base pairs, and theoretical genome coverage in culture-positive versus culture-negative samples (Table 2).

Concordance was summarized between culture of target BRD bacteria with phenotypic AST based on MICs (Supplementary Table S3) and identification of known ARGs (Supplementary Tables S4, S5) for unenriched samples, where ARGs were detected and for all samples enriched for 10 h and 14 h (Tables 3–5). The most common types of phenotypic resistance based on CLSI breakpoints were for the macrolides gamithromycin and tulathromycin, with only one tested sample with a tetracycline resistance isolate; all three antimicrobials are used in BRD management (Brault et al., 2019) (Tables 3–5, Supplementary Table S3). Samples with sulfadimethoxine MICs >256 μg/mL were also considered in the analysis based on the prevalence (Supplementary Table S3) and the use of trimethoprim-sulfamethoxazole in the treatment of BRD (Brault et al., 2019).

Detection of ARGs in unenriched samples was limited to ≤3 reads in two samples (Tables 3–5). In one sample, where the unenriched theoretical coverage of M. haemolytica was 2.8 (Supplementary Table S2A, sample 2045Bi2-023), and AST phenotypes included macrolide resistance and sulfadimethoxine MICs >256 μg/mL, 2 sul2, 2 mphE, 1 msrE, and 1 tet(34) genes were detected in the unenriched sequence data (Table 3, Supplementary Tables S4, S5A). In the second sample (2045Bi2-003), which was culture-positive for H. somni (Supplementary Table S2C) with no detected phenotypic resistance (Table 5), a single tet(H) gene was identified in the unenriched sequence data (Table 5, Supplementary Tables S4, S5C).

Resistance genes were detected more frequently in enriched samples than in unenriched samples, particularly where phenotypic resistance was most prevalent (Tables 3–5). The differences between unenriched and enriched samples were significant for macrolide and sulfonamide resistance genes within M. haemolytica reads (Table 3) and for tetracycline resistance genes within P. multocida reads (Table 4). Furthermore, ARGs were more likely to be detected (> 1 ARG) regardless of enrichment status with increased total number of base pairs for M. haemolytica (p < 0.001), P. multocida (p < 0.001), or H. somni (p = 0.03).

Macrolide ARGs were identified within at least two M. haemolytica reads (Table 3, Supplementary Tables S4, S5A) at 10h (11/20 samples) and 14 h (8/20 samples). Agreement between detection of at least two M. haemolytica reads with macrolide ARGs with phenotypic AMR for gamithromycin and tulathromycin was substantial (κ = 0.68) for 10 h and almost perfect (κ = 0.83) for 14 h (Table 3). The number of M. haemolytica reads with macrolide ARGs was also higher at 10h (p = 0.02) and 14 h (p = 0.01) in samples where M. haemolytica isolates displayed phenotypic macrolide resistance compared with samples susceptible M. haemolytica isolates.

Genes coding for macrolide resistance were identified in 463 unique M. haemolytica reads sequenced in the enriched samples by at least one of the NCBI, CARD, or MEGARes databases (Supplementary Tables S4, S5A). The 859 macrolide genes on M. haemolytica reads had a median identity of 96% (IQR, 93–97%) and a median coverage of 98% (IQR, 97–99%) (Supplemental Table S5A). The macrolide resistance-associated genes identified on M. haemolytica reads in these samples were mphE (53%, 457/859) and msrE (47%, 402/859). Both mphE and msrE were identified on 81% of the 463 reads that contained at least one of these genes. Both genes were also identified on M. haemolytica reads in all eight samples with macrolide-resistant isolates and five of the remaining eight samples with susceptible isolates.

Macrolide ARGs were also identified within at least two P. multocida reads (Table 4, Supplementary Tables S4, S5B) from samples at 10h (3/20 samples) and 14 h (5/20 samples). However, there were no samples with cultured macrolide-resistant P. multocida (Table 4, Supplementary Table S3). Macrolide ARGs (53% msrE (25/47) and 47% mphE (22/47)) were identified on 24 P. multocida reads in enriched samples by at least one database. The median gene identity was 96% (IQR, 89–99%) and the median coverage was 98% (IQR, 95–100%) (Supplementary Table S5B).

No macrolide resistance-associated genes were identified within the H. somni reads from the 20 enriched samples (Table 5, Supplementary Tables S4, S5C). None of the samples had H. somni isolates that displayed phenotypic resistance to macrolides (Table 5, Supplementary Table S3).

Macrolide ARGs (msrE, mphE) were identified in 46 B. trehalosi reads from four samples enriched for 10 h and the same four samples plus two more enriched for 14 h. Both genes had been identified in M. haemolytica reads from the same six samples.

M. haemolytica reads containing sulfonamide ARGs (Table 3, Supplementary Tables S4, S5A) were identified at least twice in samples for 10 h (12/20 samples) and 14 h (13/20 samples). Agreement for detection of at least two M. haemolytica reads containing sulfonamide ARGs (sul2 gene) and samples culture-positive for M. haemolytica with MIC >256 μg/mL for sulfadimethoxine were substantial (κ = 0.65) for 10 and 14 h (Table 3, Supplementary Table S3). The number of sul2 genes detected in M. haemolytica reads was higher for samples with M. haemolytica isolates with sulfadimethoxine for MICs >256 μg/mL than those without sulfadimethoxine in 10 h-enriched samples (p = 0.03) but not in 14-h-enriched samples (p = 0.11). The sul2 gene (n = 761) was identified in 735 M. haemolytica reads from enriched samples. The median sul2 gene identity was 95% (IQR, 92–97%) and the median coverage was 98% (IQR, 97–99%).

Sulfonamide ARGs within P. multocida reads (Table 4, Supplementary Tables S4, S5B) were identified in at least two reads for samples for 10h (3/20) and 14 h (4/20). Agreement between detection of sul2 genes and phenotypic MICs >256 μg/mL was moderate for 10h- (κ = 0.60) and 14h-enriched samples (κ = 0.48) (Table 4). There was no significant difference in the number of sul2 genes detected for samples with P. multocida isolates with MICs for sulfonamides >256 μg/mL and those without MICs (Table 4). Forty sul2 genes were identified in 28 P. multocida reads from enriched samples; median identity was 95% (IQR, 93–97%) and median coverage was 98% (IQR, 97–99%).

No sulfonamide ARGs were detected within H. somni sequence data from samples enriched for either 10h or 14 h (Table 5, Supplementary Tables S4, S5C). However, the sulfadimethoxine MIC for all four H. somni isolates was >256 μg/mL.

A sulfonamide ARG (sul2) was identified within B. trehalosi in 56 reads from four samples enriched for 10 h and the same four samples plus five additional samples enriched for 14 h.

Tetracycline ARGs were detected within at least two M. haemolytica reads (Table 3, Supplementary Tables S4, S5A) in 10h (11/20) and 14h-enriched (14/20) samples; however, there was no detected phenotypic resistance in culture-positive isolates (Table 3, Supplementary Table S3). In all enriched samples, 277 tetracycline ARGs were detected in 275 M. haemolytica reads. The primary tetracycline ARGs included tet(H) (60%, 167/277) and tet(34) (39%, 109/277). The median identity for the tet(H) genes was 95% (IQR, 93–97%) and median coverage was 98% (IQR, 97–99%) but median identity for the tet(34) genes was only 62% (IQR, 61–62%) and median coverage was 77% (IQR, 67–84%). The tet(34) genes were not detected by the Abricate option and MEGARes databases likely due to the higher cutoff of 80% identity. The tet(H) genes were detected in the same six samples for 10h and 14 h.

Tetracycline ARGs were detected within at least two P. multocida reads (Supplementary Tables S4, S5B) for 10h (4/20) and 14 h (6/20) (Table 4). The tet(H) gene was detected in at least two P. multocida reads in four samples for 10 h and five samples for 14 h. One sample contained two tet(34) genes with low identity scores (62%) as reported for the same gene from M. haemolytica. No P. multocida isolates were phenotypically resistant to tetracycline (Table 4, Supplementary Table S3).

The tet(H) genes were detected in at least two H. somni reads for two samples for 10h and 14 h. The tet(H) gene was identified in one of four H. somni culture-positive samples with tetracycline-resistant isolates (Table 5, Supplementary Tables S3, S4, S5C).

The tetracycline resistance gene (tet(H)) was identified within B. trehalosi in eight reads from one sample for 10 h and the same sample enriched for 14 h.

The only other ARGs reported at least twice in enriched samples conferred aminoglycoside resistance. The most common was APH(3′). There were 282 APH(3′)-Ia genes identified on 278 M. haemolytica reads from 14 samples and 13 genes on 12 P. multocida reads from 4 samples (Supplementary Tables S4A,B, S5A,B). There were also 232 APH(3′)-Ib genes identified on 230 M. haemolytica reads from 14 samples, and 20 genes were identified on 18 P. multocida reads from eight samples.

Another identified aminoglycoside resistance gene was APH(6)-Id, with 245 APH(6)-Id genes identified on 240 M. haemolytica reads from 14 samples and 13 genes on 12 P. multocida reads from 4 samples (Supplementary Tables S4A,B, S5A,B).

No aminoglycoside ARGs were identified in H. somni reads (Supplementary Tables S4A,B, S5C) or B. trehalosi reads.

Multiple ARGs were detected on 739 long individual reads for 31 metagenomic analyses of 19 unique samples (Supplementary Table S6). In total, 2 of the 739 reads with multiple ARGs were from unenriched samples.

The median number of unique ARGs on each read was three (5th percentile = two, 95th percentile = four). The median length of individual reads on which multiple ARGs were found was 21,941 bp (5th percentile 7,971, 95th percentile 53,678 bp). Most reads with multiple ARGs were M. haemolytica (n = 694), with smaller numbers of P. multocida (n = 42) and only one H. somni read. The most common pattern of multiple ARGs per read was msrE, mphE, and sul2 in 349 reads followed by APH(3″)-Ib, APH(6)-Id, APH(3′)-Ia, and sul2 in 188 reads. Both macrolide ARGs and sulfonamide ARGs were detected in 45% of the 774 M. haemolytica reads with at least one of these genes and in at least one M. haemolytica read from 78% (14/18) of samples with at least one of these genes.