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) via antigen capture ELISA, 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 (Additional file 4: 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 (Additional file 4: Table S2) [25]. Cattle were removed from pens if they had a clinical score of ≥ 2. Animals were classified as BRD positive if they had a rectal of temperature ≥ 40 °C and/or a clinical score of ≥ 3. Animals were treated for BRD with antimicrobials based on the feedlot protocol (Additional file 4: 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 [17], 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 [11]. 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); swabs were inserted approximately 2–3-inch into the nasal passages for sampling. 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 by passing swabs to the caudal limit of the nasopharynx at the level of the palatopharyngeal arch; 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 to stabilize the microbial community structure 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.
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 h 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 h, 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 h 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 h 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 [26].
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 [27]. Isolates were characterized as multidrug resistant (MDR) if they were not susceptible to antimicrobial(s) from ≥ 3 antimicrobial classes [28]. 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
DNA was isolated from swab samples using a QIAamp PowerFecal DNA Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. Following isolation, DNA was quantified (ng/uL) using a Qubit Flex fluorometer (ThermoFisher, Waltham, MA).
qPCR sample preparation and reaction conditions
From the extracted DNA, two 400 ng DNA aliquots were sent to Mississippi State University 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 a Qubit 4 fluorometer (Thermofisher, Waltham, MA), and the mean concentration sample DNA templates was 6.73 ± 2.00 ng/μL. A standard curve for DNA quantification was made using 8, tenfold dilutions (maximum = 1.8 ng/μL, minimum = 1.8 × 10−7 ng/μL) of DNA extracted from a pure culure of Mh confirmed by Sensititre GNID (Thermofisher, Waltham, MA). Nine replicates of each standard were made and run using a QuantStudio 3 Real-Time PCR instrument (Thermofisher, Waltham, MA) with the following reaction mixture: 7 μL of Mh DNA standard, 10 μL of PerfeCTa SYBR Green FastMix Low ROX (Quantabio, Beverly, MA), 1 μL each of forward (F) and reverse (R) primer for Mh leukotoxin D gene (lktD) (F-CTGCAACAAAGCCGATATCTTT, R-TACGACTGCTGAAACCTTGAT) [12], and molecular grade H2O to reach a final volume of 20 μL. Amplification occurred under the following conditions: 95 °C for 5 min, then 45 cycles of 95 °C for 15 s, and 60 °C for 45 s. QuantStudio Design and Analysis Software v. 1.5.1 default settings (Thermofisher, Waltham, MA) were used to determine cycles of quantification (Cq) threshold of replicates, melt-curve analysis, and other quality control checks, and results were exported to a spreadsheet for analysis using Excel for Mac Version 16.5 (Microsoft). A standard curve of Cq v. log10(ng of DNA) was created for calculation of the mass of Mh DNA. The lowest mass of DNA with SD Cq ≤ 0.5 was considered the limit of detection.
For sample plates, all reactions 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 sample DNA, 10 μL of PerfeCTa SYBR Green FastMix Low ROX (Quantabio, Beverly, MA), 1 μL each of F and R primer for Mh lktD gene [12], and molecular grade H2O to reach a final volume of 20 μL. A smaller calibration curve using five, tenfold dilutions of DNA extracted from pure growth of Mh confirmed by Sensititre GNID (Thermofisher, Waltham, MA) were included on each 96-well MicroAmp plate (4316,813, Thermofisher, Waltham, MA) to confirm reaction efficiency between 90–110%. Any plates with reaction efficiency less than 90% or greater than 110% were rerun. Also included on each plate were 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 same conditions as described above.
Cq was determined using QuantStudio Experiment Design and Analysis Software v. 1.5.1, then reviewed manually. Melt curves were used to check reaction specificity. Samples with undetermined Cq, with Cq SD greater than 0.5, with melt curves indicating non-specific binding, and/or with calculated DNA mass of less than limit of detection determined from overall standard curve, were considered to have no amplification. Mass of Mh DNA was calculated from standard curve and was logarithmically transformed for statistical analysis of geometric means. For samples with no amplification, the mass of Mh DNA was recorded as 1 × 104 ng, a non-zero number below the limit of quantification. Log10 (Mh DNA) per ng of DNA in reaction was recorded and used as outcome variable for statistical analysis.
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) [29]. The resulting pooled amplicon library was sequenced on an Illumina NovaSeq instrument using paired-end chemistry (2 × 250 bp) 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 [30]. Amplicon sequence variants (ASVs) were generated using DADA2 [31], 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 [32], 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 [33] 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 Additional file 4: Table S3. ASV counts for each sample were then normalized using cumulative sum scaling [34] and beta-diversity was analyzed using a generalized UniFrac distance matrix [35, 36]. 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 [37] and pairwiseAdonis [38] 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 [39]. 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) using the rstatix and stats packages in R [39, 40]. 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 [40].
Comparisons of log10(ng Mh DNA) per nanogram of DNA among swab types and Mh culture status were assessed using Kruskal–Wallis analysis of variance by ranks using rstatix, stats, and diplyr packages [39,40,41]. If differences were found, pairwise comparisons were tested with a Dunn test with Benjamini–Hochberg correction for multiple comparisons. Differences in qPCR amplification (Yes or No) rates among 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.