Animal experiment and sampling
Animal experiments were completed at the University of Saskatchewan following guidelines provided by the Canadian Council on Animal Care and approved by the University of Saskatchewan Animal Care Committee (Protocol #20170015). Calves recruited to the study were 5- to 6-month-old, suckling Hereford-cross females reared within the same herd (Goodale Farm, University of Saskatchewan, Saskatoon, SK Canada). Thirty (30) calves were randomly assigned to three treatment groups using Tufts Randomization plan. Experimental groups were: Suckling calves—calves remained with dams; Weaned—calves separated from dams on experimental day 0; Weaned + Transport—calves separated from dams on experimental day 0 (D0) and transported for 4.5 h before returning to the same research facility (Fig. 6). Calves in the Weaned and Weaned + Transport group were co-mingled in a drylot with access to water and brome-alfalfa hay throughout the trial. Suckling calves remained with their dams in an adjacent paddock with access to water and brome-alfalfa hay throughout the trial. On D0, blood and deep nasal pharyngeal swabs (DNS) were collected from all calves immediately prior to separation of the Weaned and Weaned + Transport calves from their dams (Fig. 6). Blood samples for serum and leukocyte isolation was collected from the jugular vein using 10 ml BD Vacutainer™ SST and EDTA blood collection tubes (Becton Dickenson, Franklin Lakes, NJ). DNS were collected using double guarded culture swabs (Jorgensen Laboratories Inc., Loveland CO). The guarded swab was inserted into the nasal cavity a distance approximately equal to the distance from the external nares to the medial canthus of the eye. The sterile swab was then extended beyond the sheath until it contacted an obstruction and the swab was rotated three times against the mucosal surface. The swab was retracted into the protective sheath, removed from the nostril, placed in a sterile tube, transported on ice to the lab, and stored at -80 °C until DNA was extracted. Blood and DNS were collected again from all calves on days 2, 4, 8, 14, and 28 following separation of the Weaned and Weaned + Transport calves from their dams (Fig. 6). Duplicate one ml aliquots of serum and two snap-frozen pellets of 10 million blood leukocytes were stored at -20 °C and -80 °C, respectively, until analyzed.
Isolation of blood leukocytes
Briefly, 12 mL of erythrocyte lysis buffer (0.17 M NH4Cl, 10 mM KHCO3, and 0.11 mM EDTA; pH 7.3) was added to 3 mL bovine blood. Cells were centrifuged at 325 g for 8 min and the supernatant discarded. Cell pellets were re-suspended in 1 mL Dulbecco’s Modified Eagle Medium (Sigma Aldrich) containing 10% fetal bovine serum (FBS) and cells counted with a Z2 Coulter Particle Count and Size Analyzer (Beckman Coulter, Brea, CA). Aliquots of 10 million blood leukocytes were pelleted at 311 g for 8 min and cell pellets snap-frozen in liquid nitrogen and stored at − 80 °C.
RNA isolation from blood leukocytes
RNA was extracted from blood leukocytes using a combined TRIzol/RNEasy Mini Kit extraction method. Frozen cell pellets were suspended in 1 mL TRIzol reagent (ThermoFisher Scientific, Waltham, MA) and 200 µL chloroform (Sigma Aldrich, St. Louis, MO) was added to each sample. Samples were shaken for 15 s and incubated at room temperature for 2–3 min before centrifuging for 15 min at 13,282 g. Following centrifugation, the aqueous phase was removed, an equal amount of 70% ethanol added, and samples were added to the silica columns provided in the RNEasy Mini Kit (Qiagen, Hilden, North Rhine-Westphalia, Germany) and processed according to manufacturer’s instructions.
cDNA synthesis and reverse transcription PCR
Synthesis of complementary DNA (cDNA) from blood leukocyte RNA template was performed following manufacturer’s instructions for the Quantitect Reverse Transcription Kit (Qiagen, Hilden, North Rhine-Westphalia, Germany). A 30-min cDNA synthesis incubation step (42 °C) was used to remove excess RNA secondary structure. Briefly, 500 ng RNA was diluted in 6 µL UltraPure DNAse/RNAse-Free Distilled Water (Invitrogen) and 1 µL of 7 × gDNA wipe-out buffer was added to remove genomic DNA. The GeneAmp 9700 PCR System (Applied Biosystems, California USA) was used to incubate this mixture for 2 min at 42 °C. Following incubation, 3 µL master mix was added to each reaction. The master mix for each reaction consisted of 2 µL 5 × Quantiscript RT buffer, 0.5 µL of primer mix, and 0.5 µL reverse transcriptase. Following addition of the master mix, each reaction was incubated for 30 min at 42 °C, followed by 3 min at 95 °C.
Reverse transcription PCR (RT-PCR) reactions were prepared with 25 ng of cDNA (5 µL of 5 ng/µL cDNA) and 10 µL of master mix. The master mix consisted of 7.5 µL 2 × PERFECTA-IQ SYBR Green Supermix (QuantaBio), 2.2 µL UltraPure DNAse/RNAse-Free Distilled Water (Invitrogen, Carlsbad, CA), and 0.3 µL of 10 µM forward and reverse primer (3 pmol; Additional file 2: Table S4). Reactions were run in Hard Shell Low-Profile 96-well semi-skirted, clear-shell, and clear-well PCR plates (BioRad, Hercules, CA). The CFX Connect Real Time System (BioRad, Hercules, CA) was used to run and quantify the real time PCR reactions. Reactions were first run at 95 °C for 2 min to activate the hot-start Taq polymerase, then for 40 cycles at 95 °C for 15 s (denature), 60 °C for 30 s (anneal), and 72 °C for 30 s (extend). Following amplification, a melt curve was applied for detection of abnormal products. The melt curve started at 65 °C, and the temperature held for 10 s before increasing by 1 °C. This pattern was repeated to a temperature of 95 °C. Results were visualized using the CFX Manager/Maestro software and corrections for primer efficiency were included in Cq value calculations.
Serum antibody titres for M. heamolytica, P. multocida, and My. bovis
Antibody capture enzyme-linked immunosorbent assays (ELISAs) were performed to quantify serum IgG antibody titers (Fig. 6) using the protocol described in Hill et al. [22]. Multiple Mycoplasma spp. are present within the bovine URT but the antibody capture ELISA used was specific for My. bovis [28]. Briefly, proteins used for antibody capture included recombinant M. haemolytica leukotoxin [29], recombinant My. bovis MilA-ab [28] and soluble bacterial lysates prepared from M. haemolytica and P. multocida. Serum titers are presented as the inverse of the final serum dilution generating an OD reading exceeding the mean + 2 SD of the background OD value from triplicate wells containing negative serum samples.
Profiling the URT microbiome using metagenomics sequencing
Total DNA was extracted from the DNS using PowerSoil DNA isolation kit (MO BIO Laboratory Inc., Carlsbad, CA). Briefly, the DNS was transferred into a PowerBeads tube containing the C1 solution and subjected to bead-beating using Mini-BeadBeater-16 (BioSpec Products, Bartlesville, OK) at 5000 rpm for 3 min. The supernatant was separated after centrifuging at 13,000 rpm for 15 min and subsequently used to isolate DNA following the manufacturer’s instructions. DNA quantity was measured using Qubit 3.0 fluorometer (Invitrogen, Carlsbad, CA) and Qubit dsDNA HS assay kit (ThermoFisher Scientific, Waltham, MA). Shotgun DNA libraries (Fig. 6) were prepared using NEB Ultra II DNA Library Prep Kit for Illumina (New England Biolabs Inc., Massachusetts, USA) and sequenced using Illumina HiSeq4000 PE100 (Illumina, California, USA) at Genome Quebec (McGill, Quebec).
Analysis of metagenomics sequencing data
Demultiplexed raw data (229.8 Gb, Additional file 2: Table S5) were first run through Trimmomatic version 0.39 [30] in paired-end mode to remove adapters, low quality sequences (Phred < 20) and short sequences (< 50 bp). Then, host contamination was removed using Bowtie 2 [31], SAMtools [32] and BEDtools [33] by aligning sequences to bovine genome (UMD 3.1). Unassembled sequences with host contamination removed (53.3 Gb) were then uploaded into the MG-RAST metagenomic analysis server [34], version 4.0, and paired-ends were joined for each sample before submitting for processing. Artificial replicates, host (bovine) DNA and low-quality (Phred score < 25) sequences were removed from the raw data, and the remaining good-quality sequences were used to assign the microbial functions using the subsystems annotation source in the SEED hierarchy and KEGG Orthology and microbial taxa using RefSeq database. A maximum cut-off e-value of 1e-10, maximum identity of 70% and maximum alignment length of 80 was used as data selection criteria for the functions and taxa abundance analyses. In addition, MEGAHIT v1.1.1 [35] was used to assemble raw sequences with a minimum contig length of 200 bp and K-mer size 119. Assembled sequences were then used to assign taxonomic composition using Kraken2 database [36,37,38].
Estimation of total bacteria M. haemolytica and P. multocida densities
The density of total bacteria, M. haemolytica and P. multocida were estimated using quantitative real-time PCR (qPCR) and bacterial primers (Additional file 2: Table S6). Total DNA extracted from DNS was diluted to 50 ng/µL and 1 µg template was used to perform qPCR with SYBR Green chemistry (Fast SYBR® Green Master Mix; Applied Biosystems, Foster City, CA) and StepOnePlus™ real-time PCR system (Applied Biosystems, Foster City, CA). The standard curve of total bacteria was constructed using purified PCR products amplified with 27F and 1492R primer pair, while standard curves for M. haemolytica and P. multocida were constructed using genomic DNA extracted from pure cultures of each bacterial species. Bacterial density (copy number of the 16S rRNA gene per DNS for total bacteria and P. multocida and copy number of leukotoxin (lkt) gene per DNS for M. haemolytica) was calculated using the following equation: (quantity mean × DNA concentration × DNA elution volume)/DNA amount used in qPCR reaction.
Statistical analysis
The ADR gene expression data and antibody ELISA data were analyzed using a repeated measure model with sampling point as the repeated measure and the generalized least square function using autoregressive of order 1 (AR1) covariance structure. Significant time and treatment effects were observed and a post-hoc test for multiple comparisons of factors was performed using Tukey’s multiple comparison test to determine if there were significant treatment effects.
Taxonomic (at genus level) and functional (at KEGG Orthology level 2) compositions all metagenomes were first analyzed using principal component analysis (PCA) and Bray–Curtis dissimilarity matrix to understand the effect of stressors and sampling time point on URT microbiome. Analysis of similarities (ANOSIM) was used to test statistical significances of the PCA-based visualization. Then, the non-parametric Kruskal–Wallis test one way ANOVA by rank was performed to test the effect of stressors and sampling time point on the relative abundance of potentially pathogenic bacterial groups (Mannheimia, Pasteurella, Haemophilus, Histophilus, Moraxella, Mycoplasma) and bacteriophage (Microvirus, P2-like viruses). The same analysis was performed after stratifying data by stressor type to compare time points to test interaction effect between type of stressor and time point. A post-hoc test for multiple comparisons of factors was performed using pairwise Wilcox (Mann–Whitney U-tests) test and p values were adjusted using Benjamini and Hochberg method [31]. Data were presented as medians with 95% confidence intervals (CIs) and statistical differences were declared at p-value adjusted < 0.05.
A logistic regression analysis was performed to explore the relationship between the colonization of potential pathogens and type of stressor. High abundance of potential pathogens in the logistic regression model was defined as above (Yes) and below (NO) the median relative abundance and used to calculate the odds ratio (OR) of colonization (R package “questionr”). A stratified logistic regression model was fitted to identified temporal effect on the likelihood of colonization of microbial taxa within in each treatment group. The effect of stress on the likelihood of colonization of microbial taxa was declared in relative to suckling group after adjusting the model for sampling time. While the effect of sampling time was declared relative to D0 and D2. Associations between host immune responses to M. haemolytica and P. multocida and the abundance of genera Mannheimia and Pasteurella in the URT were explored using two different approaches. First, a negative binomial regression model was fitted for three groups together and within each stressor type. The colonization of Mannheimia and Pasteurella was defined as above (Yes) and below (NO) the median relative abundance. A mediation analysis (R package “mediation”) was then performed to test if changes in host immune responses were mediated by the respective bacterial group colonized in the URT. Bacterial densities estimated through qPCR were first normalized by log10 transformation and day 0 data were then analyzed using a one-way ANOVA to test effect of calf group on initial bacterial densities. Post-weaning data (days 2, 4, 8, 14 and 28) were then analyzed using a repeated measure model with sampling point as the repeated measure and the generalized least square function using autoregressive of order 1 (AR1) covariance structure, which was selected as the best fit by the Bayesian information criterion. All data were analyzed using R package (version 4.0.0).