Animal populations and sample collection
A total of eight, six to eight month old, single-source, Charolais feedlot calves (mean body weight 348.47 ± 14.59 kg), were enrolled in our study in April 2016, approximately 60 days after arrival at the South Farms-Beef Cattle and Sheep Field Laboratory, University of Illinois at Urbana- Champaign- USA. All calves were clinically healthy, with no history of receiving any antimicrobial drugs prior to, or after arrival at the feedlot. The calves were vaccinated immediately after arrival at the feedlot, with a modified live virus vaccine against IBR, BVDV (Types I and II), BPIV-3 and BRSV (5 ml IM; Bovi-Shield Gold FP5 L5 HB Cattle Vaccine, Zoetis Animal Health), and dewormed with a topical anthelminthic (Noromectin Pour-On Solution, Norbrook® Inc. USA). The calves were fed a mixed ration (20% silage, 20% modified wet distillers grains with soluble, 10% dry supplement, and 50% high moisture corn), and given free access to water. The use of animals for this study was approved by the University of Illinois Institutional Animal Care and Use Committee (IACUC Protocol: #15064) and all of the experimental protocols were performed in accordance with relevant guidelines and regulations set by IACUC.
Physical exams, including rectal temperature, respiratory rate and pulse rate, were performed with the calves standing in the hydraulic chute. A clinical lung score was recorded using the microphone of an automated Whisper stethoscope (Whisper®, Geissler Corp, Plymouth, MN). Briefly, the stethoscope was placed over the 5th intercostal space of the right thoracic wall, approximately 10 cm above the elbow, at a site known to encompass the apical lung lobes. Recorded lung sounds were then automatically transmitted wirelessly to a computer located within 2 m of the stethoscope, and analyzed using software rendering a 5-point lung score scales [19]. The program software is designed to remove heart sounds and potential interference from the environment, and classifies acoustic patterns in to lung scores ranging from 1 to 5 (1 = normal, 2 = mild acute, 3 = moderate acute, 4 = severe acute, and 5 = chronic).
Feedlot calves with a rectal temperature < 39.4 °C, respiratory rate < 50 breaths/min, pulse rate < 120 beats/min and lung score ≤ 2 were defined as healthy, and selected for sampling.
Following physical examination, a single, deep NPS was collected from each calf, using a 33-in.-long double guarded PVC culture swab (Kalayjian Industries, Inc. U.S.A.) according to published techniques [20]. Briefly, the calf’s nostrils were cleaned with a disposable wipe before collection. The NPS was carefully inserted into the ventral meatus of the nose, and advanced approximately 2/3 of the dorsal head length. Once in place, the end of the swab was exposed to a length of approximately 1–2 in., by withdrawing the outer sleeve, and the sampling unit was then firmly rotated through 360° against the pharyngeal wall, for 20–30 s. The cotton tipped swab was then retracted back in to the outer sheath, and the whole swab was removed gently from the animal’s nose. Each cotton swab was then broken off into a sterile 2 ml cryo-tube, transported on dry ice to the laboratory, and stored at −20 °C pending further processing.
Following NPS sampling, the calves were sedated with Xylazine hydrochloride (0.1 mg/kg IM) (Rompun®, Bayer health care LLC, Kansas) and the nostrils were cleaned with dry gauze. BAL was performed with the sedated calves in a standing position, using a sterilized, flexible BAL tube (BAL300 – 300 cm length, 10 mm outside diameter, 2.5 mm internal diameter, balloon volume 10 cm3, MILA International, Inc. U.S.A.) [21]. The BAL catheter was introduced into the nostril, directed into the ventral meatus, and then advanced until it encountered resistance in the caudal pharynx. At this point, the calf’s head and neck were held straight, in maximal extension, to allow the catheter to pass into the trachea during the inspiratory phase of the respiratory cycle. On reaching a wedged position in the lower respiratory airway, the catheter was secured by inflating the balloon cuff with 5 cm3 of air. For sampling, 120 mls of sterile saline was infused through the catheter lumen, using 60 ml syringes attached to a stopcock and catheter tipped adapter. Immediately after the 120 ml infusion, negative pressure was applied to aspirate airway fluid, which was immediately placed into a sterile 50 ml specimen tube. The BAL samples were stored on ice until processing, approximately 2–4 h after collection. The collected BAL was centrifuged at 14,000×g for 10 min. The pellets were re-suspended in 500 μL sterile PBS and stored at −20 °C pending further analysis.
Genomic DNA extraction
Genomic DNA extraction was performed from each NPS and BAL sample using the PowerFecal® DNA isolation Kit (MO BIO Laboratories, Inc., Carlsbad, CA, USA), according to the manufacturer’s instructions. Briefly, each sample was transferred to a dry bead tube with 60 μl of Solution C1 and 750 μl of Bead Solution, heated at 65 °C for 10 min, and settled in Bullet Blender 24 Gold tube holder (Next Advance, Inc., Averill Park, NY, USA). The tubes were vortexed at maximum speed for 10 min to achieve microbial cell disruption. The PowerFecal® DNA isolation Kit protocol was used to complete the extraction, according to manufacturer instructions. Purified DNA was eluted into 50 μl of Solution C6 rather than 100 μl to increase the DNA concentration. Total DNA concentration was quantified using a Nanodrop™ spectrophotometer (NanoDrop Technologies, Rockland, DE, USA) at wavelengths of 260 and 280 nm, and the integrity was confirmed by agarose gel electrophoresis.
Fluidigm access Array amplification of the V3-V4 hypervariable region of 16S rRNA genes and Illumina sequencing
Genomic DNA was subject to Fluidigm Access Array Amplification (Fluidigm Corporation, South San Francisco, CA, USA). Prior to amplification, all DNA samples were measured on a Qubit™ fluorometer (Life technologies, Grand Island, NY, USA) using the High Sensitivity DNA Kit (Agilent Technologies, Santa Clara, CA, USA) [22].
Briefly, the primer sequences F357-for (CCTACGGGAGGCAGCAG) and R806-rev (GGACTACNVGGGTWTCTAAT) were used to amplify the V3-V4 hypervariable region of the 16 s rRNA gene. PCR reactions were performed on a Fluidigm Biomark HD™ PCR machine (Fluidigm Corporation, South San Francisco, CA, USA) using the default Access Array cycling program without imaging.
Harvested product was quantified on a Qubit™ fluorometer (Life technologies, Grand Island, NY, USA) and stored at −20 °C. All samples were run on a Fragment Analyzer™ (Advanced Analytics, Ames, IA, USA), and the amplicon regions were quantified. PCR products were then size selected on a 2% agarose E-gel™ (Life technologies, Grand Island, NY, USA) and extracted from the isolated gel slice with the QIAquick Gel extraction kit (Qiagen, Valencia, CA, USA). Cleaned, size-selected product was examined on an Agilent Bioanalyzer™ to confirm appropriate profile, and for the determination of average size. The final pooled Fluidigm libraries were transferred to the DNA Services lab at the W. M. Keck Center for Comparative and Functional Genomics (University of Illinois at Urbana-Champaign, Urbana, IL, USA) for Illumina sequencing. The Illumina MiSeq™ platform (Illumina, San Diego, CA, USA) was used to sequence the V3- V4 region of the 16S rRNA gene according to the Illumina instructions.
Sequence data processing and statistical analysis
The raw sequence data were preprocessed and analyzed using the open-source software package, Quantitative Insights Into Microbial Ecology (QIIME®) version 1.9 (http://qiime.org/) [23]. Sequences were filtered for quality the default parameters of the split_libraries.py command; minimum sequence length equal 200, maximum sequence length equal 1000, a Phred score of less than 25, maximum number of ambiguous bases equal 6 and homopolymer runs of >6 bp [24]. Chimeric sequences were detected and removed using UCHIME [25]. The remaining sequences were clustered into operational taxonomic units (OTUs) using open reference OTU selection protocols (97% identity cutoff) with the UCLUST algorithm [26], and assigned a taxonomic classification against the Greengenes® database [27].
The core microbiota, those microbes shared among all sampled calves, was identified at the genus level. For subsequent alpha and beta diversity analysis, the OTU table was randomly subsampled and rarefied to 2400 sequences per sample. Alpha diversity (an estimate of bacterial community richness in a sample) was calculated within QIIME® using the Chao1 (an estimate of species richness), observed species (the total number of microbial species present in a community), and phylogenetic diversity (PD whole tree) (an estimate of the biodiversity which incorporates phylogenetic difference between species). Beta-diversity (an estimate of bacterial communities expression of diversity between different sites) was calculated using weighted and unweighted UniFrac distance [28] and clustering of the samples based on distance was visualized on principal coordinate analysis (PCoA) plots using EMPeror® [29].
The OTU relative abundance values were analyzed using the linear discriminant analysis (LDA) effect size (LEfSe) algorithm, to identify OTUs that display significant differences between the two sites [30]. Additionally, a cladogram was produced using the online LEfSe tool. The algorithm first used the non-parametric factorial Kruskal-Wallis test to detect taxa with significantly different abundance, followed by pairwise Wilcoxon test to detect biological consistency between NPS and BAL samples, and then used LDA to estimate the effect size of each differentially abundant feature. Statistical analyses were performed using JMP® 12.12 (SAS Institute Inc., North Carolina, USA).
For characterization of NPS versus BAL microbiota, analysis of bacterial abundance between the NPS and BAL samples in healthy calves was performed using nonparametric Wilcoxon tests. Alpha diversity metrics (Chao1, observed species and PD whole tree) were also compared between the two groups using a non-parametric Wilcoxon tests in JMP® 12.12 software. The weighted and unweighted UniFrac distances for NPS and BAL samples were also compared using nonparametric ANOSIM test (analysis of similarities) with 999 Monte Carlo permutations.
To further explore the possible relationships between upper and lower respiratory tract microbial communities, multivariate analysis (linear correlation matrixes) was performed using JMP® 12.12 (SAS Institute Inc., North Carolina, USA). In view of the large number of comparisons, an additional, more conservative approach was adopted using multiple linear forward regression analysis. Forward linear regression analysis was performed using the SPSS version 22 (IBM Corp, Version 22.0, Armonk, NY, 2013) statistical package. Differences with a P value ≤0.05 were considered significant for all analyses.
Fastq data obtained in the current study were uploaded to the sequence read archive (SRA) on the National Center for Biotechnology Information (NCBI) website (http://www.ncbi.nlm.nih.gov/sra), to make the files available for a public database with bio-project accession number PRJNA323521.