Study population and sampling
Two groups each consisting of 60 beef-type steers and bulls were purchased from a livestock auction market located in central Texas, shipped to the West Texas A&M University Research Feedlot on May 14 and May 21, 2020, where they were enrolled in this study (n=120). Upon arrival at the feedlot, cattle received an ear tag with an individual identification number and were processed following standard practices of many feedlots. Briefly, tildipirosin (Zuprevo, Intervet Inc., Summit, NJ), a long-acting macrolide, was administered to every animal at 4 mg/kg subcutaneously for BRD metaphylaxis. Animals were vaccinated against clostridial (Calvary 9, Merck Animal Health, Omaha, NE) and respiratory bacterial pathogens (Once PMH, Merck Animal Health, Omaha, NE), given a zeranol growth implant (Ralgro, Merck Animal Health, Summit, NJ) and given anthelminthic therapy with albendazole (Valbazen, Zoetis, Kalamazoo, MI) and ivermectin with corsulon (Ivermectin Plus, Durvet, Inc., Blue Springs, MO). Animals were also tested to identify any animals persistently infected with bovine viral diarrhea virus (BVD-PI), and any BVD-PI animals were removed from the study. Bulls were castrated and given meloxicam at 1.1 mg/kg orally the day following metaphylaxis, vaccine, and anthelminthic administration (Table S1).
Pens were monitored daily by trained feedlot personnel to identify animals with BRD, and animals were assigned a BRD clinical score of 0-4 based on visual signs of disease (Perino & Apley, 1998; Table S2). Cattle were removed from pens if they had a clinical score of \(\ge\) 2. Animals were classified as BRD positive if they had a rectal of temperature \(\ge\) 40°C and/or a clinical score of \(\ge\) 3. Animals were treated for BRD with antimicrobials based on the feedlot protocol (Table S2). The animals were on feed for 213 and 255 days for group 1 and group 2, respectively.
On day 14 after arrival, when a high prevalence of M. haemolytica shedding was expected (Woolums, 2018), cattle were processed through a chute, where they were weighed and restrained for sampling. Six different nasal and nasopharyngeal samples (three from the left and three from the right) were obtained as previously described (Godinho et al., 2007). Briefly, the external nares were cleaned with a paper towel to remove superficial secretions and dirt, and both internal nasal passages were then swabbed with the 6-inch (15.2 cm) rayon fiber nasal swabs (NS) (SP130D, Starplex Scientific Corporation, St. Louis, MO). After collecting nasal swabs, the 16-inch (40.6 cm) rayon fiber proctology swab (PS) (816-100, Puritan, Guilford, ME) or the 29.5 in (74.9 cm) cotton fiber deep-guarded swab (DG) (E9-5200, Continental Plastic, Delavan, WI) were used to sample the left and right nasal and nasopharyngeal passages; the order of collection of the proctology and deep-guarded swabs was randomized. All swabs collected via the left nostril were placed in modified Amies transport media (Starplex Scientific Corporation, St. Louis, MO) and used for aerobic bacterial culture, identification, and antimicrobial susceptibility testing. All swabs collected via the right nostril were placed in 100% ethanol and were used for DNA extraction and subsequent analyses with 16S rRNA gene sequencing and qPCR. All samples were kept on ice and transported to the laboratory for processing immediately after collection.
The unique animal ID was incorrectly recorded for two enrolled animals, which prevented extraction of corresponding data regarding animal weight and health records. Additionally, three swab samples intended for DNA analyses (two deep-guarded swab sample and one proctology swab sample), and one deep-guarded swab intended for culture were damaged during transport to the laboratory and could not be analyzed. DNA extraction from one deep-guarded swab sample failed, as well. These data are therefore missing from the results.
3.1 Culture, microbial identification, and susceptibility testing
Swabs collected in modified Amies media were directly streaked onto one quadrant of a plate of tryptic soy agar (TSA) with 5% sheep blood (Remel, Lenexa, KS), and sterile disposable loops (Remel, Lenexa, KS) were used to streak the rest of the plate for bacterial isolation. Plates were incubated at 37°C with 5% CO2. At 24 and 48 hours of incubation, plates were monitored for growth consistent with Mh (2-3 mm, round, raised, light-grey, smooth, shiny colonies with faint β hemolysis). If colonies consistent with such growth were present, catalase, oxidase, and indole tests were performed. If preliminary biochemical tests were consistent with Mh (catalase-positive, oxidase-positive, and indole-negative), a single colony was randomly selected by choosing the Mh-like colony closest to a mark made at a random position on the bottom of the media plate and subcultured onto a new blood agar plate and returned to the incubator at the above conditions. After 24 hours, subcultures were monitored for colony phenotype and biochemical tests consistent with Mh. If present, 5-7 colonies were randomly selected with a sterile disposable loop and suspended into 1.5 mL of Brain Heart Infusion broth (B-D, Franklin Lakes, NJ) and 30 % glycerol (Thermofisher, Waltham, MA). The same loop was then used to streak one half of another blood agar plate which was then incubated as described above for 24 hours then shipped on ice to University of Nebraska-Lincoln Veterinary Diagnostic Center (UNL-VDC) to confirm identity and for antimicrobial susceptibility testing. Primary plates with no suspected Mh growth at 48 hours were considered negative for M. haemolytica.
At UNL-VDC, a single colony from the shipped plate was subcultured overnight on blood agar to ensure pure growth which was then used to confirm Mh identification and antimicrobial susceptibility testing. Matrix assisted laser desorption-ionization time-of-flight mass spectroscopy (MALDI-TOF) was used to confirm Mh identity as well as MALDI-TOF biomarker based genotyping of Mh isolates (Loy & Clawson, 2017).
Antimicrobial susceptibility testing was performed at UNL-VDC using semi-automated broth microdilution via the Sensititre system (ThermoFisher, Waltham, MA) and the bovine/porcine panel containing gamithromycin and tildipirosin (BOPO7F Vet AST Plate, ThermoFisher, Waltham, MA). Results were interpreted according to breakpoints for Mh in BRD from the Clinical and Laboratory Standards Institute (CLSI, 2018). Isolates were characterized as multidrug resistant (MDR) if they were not susceptible to antimicrobial(s) from ≥ 3 antimicrobial classes (Sweeney et al., 2018). Because the concentration range for ampicillin on the BOPO7 plate does not include CLSI breakpoints, only minimum inhibitory concentration (MIC) was recorded, and ampicillin resistance classification was not included in determination of isolates as MDR.
DNA extraction
Metagenomic DNA was isolated from swab samples using a QIAamp PowerFecal DNA Kit (Qiagen, Hilden, Germany). Following isolation, DNA was quantified (ng/uL) using a Qubit Flex fluorometer (ThermoFisher, Waltham, MA).
3.2 qPCR Sample Preparation and Reaction Conditions
From the DNA extracted DNA, two 400 ng DNA aliquots were sent to Mississippi State for qPCR. Samples from one aliquot were diluted in Low-Tris TE buffer to an estimated final concentration of 8 ng/µL. Final concentrations were measured on Qubit 4 fluorometer, and the mean DNA concentration was 6.73 ± 2.00 ng/µL.
All reactions, including samples and controls, were run in triplicate using a QuantStudio 3 Real-Time PCR instrument (Thermofisher, Waltham, MA) and the following reaction mixture: 40 ng (mean=41.5 ng, SD=4.7 ng) of metagenomic sample DNA, 10 µL of PerfeCTa SYBR Green FastMix Low ROX (Quantabio, Beverly, MA), 1 µL each of F and R primer for Mh leukotoxin D gene (lktD) (F- CTGCAACAAAGCCGATATCTTT, R- TACGACTGCTGAAACCTTGAT) (Loy et al., 2018), and molecular grade H2O to reach a final volume of 20 µL. Five positive controls of triplicate 10-fold dilutions DNA extracted from pure growth of Mh confirmed by Sensititre GNID (Thermofisher, Waltham, MA) were included on each 96-well MicroAmp plate (4316813, Thermofisher, Waltham, MA). Also included on each plate, negative controls consisting of reaction mixture of molecular grade H2O in place of template DNA. Additionally, controls with no primer added, and no master mix controls added were included. Amplification occurred under the following conditions: 95°C for 5 minutes, then 45 cycles of 95°C for 15 seconds, and 60°C for 45 seconds.
Cycle threshold (Ct) was determined using QuantStudio Experiment Design and Analysis Software, then reviewed manually. Melt curves were used to check reaction specificity. Samples containing triplicates with Ct values differing by more than two were removed, and the remaining technical duplicates were used to determine mean Ct. In order to maximize agreement between qPCR and bacterial culture, culture was used as a gold standard, and qPCR sensitivities, specificities, and accuracies were calculated at a range of Ct (25-40; Additional File 1). Samples with Ct\(\le\)30 were considered positive for Mh. For quantitative statistical analyses, samples with no amplification were considered to have Ct=50, as done previously (Loy et al., 2018).
3.3 16S rRNA library preparation, and sequencing
Preparation of libraries for sequencing of the V3-V4 region of 16S rRNA was conducted as previously described (Illumina, 2013). The V3-V4 region of the 16S rRNA gene was amplified using the 341F (5’ – CCTACGGGNGGCWGCAG – 3’) and 805R (5’ – GACTACHVGGGTATCTAATCC – 3’) primer pair (Integrated DNA Technologies, Inc, Coralville, IA) and sequencing libraries were prepared using the Nextera IDT kit (Illumina, San Diego, CA). The resulting pooled amplicon library was sequenced on an Illumina NovaSeq instrument using paired-end chemistry (2 x 250bp) at the University of Colorado Anschutz Medical Campus’ Genomics and Microarray Core.
Bioinformatics and statistics
Demultiplexed paired-end reads generated from 16S rRNA gene sequencing were imported in QIIME2 version 2020.11 (Bolyen et al., 2019). Amplicon sequence variants (ASVs) were generated using DADA2 (Callahan et al., 2016), which was also used to filter reads for quality, remove chimeric sequences, and merge overlapping paired-end reads. Forward and reverse reads were truncated at 248 bp and 250 bp, respectively. Taxonomy was assigned using a Naïve Bayes classifier trained on the Greengenes version 13_8 99% OTUs database (DeSantis et al., 2006), where sequences had been trimmed to include only the base pairs from the V3-V4 region bound by the 341F/805R primer pair. Reads mapping to chloroplast and mitochondrial sequences were filtered from the representative sequences and ASV table using the ‘filter-seqs’ and ‘filter-table’ functions, and a midpoint-rooted phylogenetic tree was generated using the ‘q2-phylogeny’ pipeline with default settings, which was used to calculate phylogeny-based diversity metrics. Data and metadata were then imported into phyloseq (McMurdie & Holmes, 2013) using the ‘import_biom’ and ‘import_qiime_sample_data’ functions and merged into a phyloseq object. The proportion of reads mapped to each taxonomic rank can be found in Table S2. ASV counts for each sample were then normalized using cumulative sum scaling (Paulson et al., 2013) and beta-diversity was analyzed using generalized UniFrac distances (Lozupone et al., 2011, Chen et al., 2012). From these distances, principal coordinates analysis (PCoA) was performed and plotted, and a permutational multivariate analysis of variance (PERMANOVA) was used to test for significant differences in community structure using the vegan (Oksanen et al., 2019) and pairwiseAdonis (Arbizu, 2017) packages. To ensure significant differences were not the result of unequal dispersion of variability between groups, permutational analyses of dispersion (PERMDISP) were conducted for all significant PERMANOVA outcomes using the vegan package. Further, the relative abundances of ASVs within each sample were calculated and plotted using phyloseq. Differences in relative abundance were tested using a pairwise Wilcoxon rank-sum test with a Benjamini-Hochberg correction for multiple comparisons in R version 3.6.0.
Summary statistics of arrival weight, number of animals treated for BRD overall and number treated at time of sampling, and days on feed (DOF) until their first BRD treatment were calculated using R version 4.0.3 (R Core Team, 2020). Comparisons between the two sampling groups were made using Wilcoxon rank-sum test for continuous outcome variables (arrival weight and DOF until first treatment) and Chi-square test for binary response variables (treatment for BRD during feeding period and treatment for BRD at the time of sampling) in the rstatix package in R (Kassambra, 2021). Cochran’s Q test was used to compare isolation of Mh by swab type using SAS software v 9.4 (SAS Institute, Cary, NC). If differences were found using Cochran’s Q test, pairwise comparisons using McNemar’s Chi-square test were performed with the rstatix package.
Comparisons of Ct between swab types and Mh culture status were assessed using Kruskal-Wallis analysis of variance by ranks using rstatix (Kassambra, 2021). If differences in Ct were found, pairwise comparisons were tested with a Dunn test with Benjamini-Hochberg correction for multiple comparisons in the rstatix package (Kassambra, 2021). Differences in qPCR positive (Ct \(\le\) 30) and negative (Ct \(>\)30) rates between swab types were tested using Cochran’s Q test in SAS software v 9.4 (SAS Institute, Cary, NC), with post hoc comparisons tested with pairwise McNemar’s Chi-square in rstatix.