Effects of treatment with enrofloxacin or tulathromycin on fecal microbiota composition and genetic function of dairy calves

The increasing concerns with antimicrobial resistance highlights the need for studies evaluating the impacts of antimicrobial use in livestock on antimicrobial resistance using new sequencing technologies. Through shotgun sequencing, we investigated the changes in the fecal microbiome composition and function, with a focus on functions related to antimicrobial resistance, of dairy calves. Heifers 2 to 3 weeks old, which were not treated with antibiotics by the farm before enrollment, were randomly allocated to one of three study groups: control (no treatment), a single treatment of enrofloxacin, or a single treatment of tulathromycin. Fecal samples were collected at days 4, 14, 56 and 112 days after enrollment, and DNA extraction and sequencing was conducted. The effect of antibiotic treatment on each taxon and genetic functional level by time (including Day 0 as a covariate) revealed few changes in the microbiota. At the genus level, enrofloxacin group had higher relative abundance of Blautia, Coprococcus and Desulfovibrio and lower abundance of Bacteroides when compared to other study groups. The SEED database was used for genetic functional analyses, which showed that calves in the enrofloxacin group started with a higher relative abundance of “Resistance to antibiotics and toxic compounds” function on Day 0, however an increase in antibiotic resistance genes after treatment with enrofloxacin was not observed. “Resistance to Fluoroquinolones” and “Erythromycin resistance”, of relevance given the study groups, were not statistically different in relative abundance between study groups. “Resistance to fluoroquinolones” increased during the study period regardless of study group. Despite small differences over the first weeks between study groups, at Day 112 the microbiota composition and genetic functional profile was similar among all study groups. In our study, enrofloxacin or tulathromycin had minimal impacts on the microbial composition and genetic functional microbiota of calves over the study period.


Introduction
There is urgent need for the judicious use of antimicrobial drugs to extend the effectiveness of currently available therapies [1,2]. Antimicrobial resistance (AMR) is a natural phenomenon, and resistance genes are present even in places that were never inhabited by humans [3,4] and in wild animals [5]. However, the use of antibiotics in human and animal medicine, as well as in plant agriculture and animal production systems, has resulted in an acceleration in the selection and spread of antimicrobial resistance [1,6]. Antimicrobial resistant enteric bacteria can be transmitted from livestock to humans through the fecal-oral route, or contamination in the food chain and environment [1,7,8]. Together with prevention of disease and infectious agent transmission, the prudent use of antibiotics is extremely important to preserve treatment effectiveness and decrease the dissemination of resistant bacteria.
Antimicrobials may be used to control and prevent the spread of the disease in production animals at high risk of developing a bacterial infectious disease. This practice is referred to as metaphylaxis. In the United States, several drugs are approved for this use in cattle, including in dairy calves, at risk of bovine respiratory disease (BRD). BRD is known to be a common cause of disease in dairy calves at 2-3 weeks of age [9,10]. In 2014 in the United States 12% of pre-weaned calves were affected with respiratory disease and almost 95% of those were treated with antimicrobial drugs [11]. Macrolides and florfenicol were the drugs of choice on 18.2 and 15.1 percent of the farms to treat BRD, followed by penicillin (8.1%) and fluoroquinolones (6.6%) [11]. In our study, we focused on enrofloxacin and tulathromycin, a fluoroquinolone and a macrolide, respectively. These are injectable, single dose antimicrobials, labeled for treatment and control of BRD.
Fluoroquinolones are broad-spectrum antimicrobial drugs that target essential bacterial enzymes DNA gyrase and DNA topoisomerase IV. Mutations in these bacterial enzymes can confer resistance, as can plasmids carrying resistance genes that protect the cells from the bactericidal effects of quinolones [12,13]. Fluoroquinolones were developed about 40 years ago and their use increased worldwide in the 1990s for treatment of Gram-negative infections in humans [14,15]. Examples of FDA-approved fluoroquinolones for use in human medicine are ciprofloxacin and levofloxacin. In 1988, the FDA approved the veterinary use of a fluoroquinolone, enrofloxacin, for dogs and cats [16]. Later, in 1998, it was approved for treating bovine respiratory disease in cattle. The use of enrofloxacin in poultry was banned in 2005 and in 2008 its use in female dairy cattle was restricted to nonlactating animals up to 20 months old. Extra-label use of enrofloxacin is strictly prohibited. Enrofloxacin has also been approved for treatment and control of swine respiratory disease since 2008. [16].
Macrolides are mainly bacteriostatic. They inhibit bacterial protein synthesis and the spectrum of action of drugs within this class varies. In human medicine, erythromycin has been available since 1952, and other current drugs in the same class are azithromycin, clarithromycin and telithromycin [17]. In veterinary medicine, erythromycin and tylosin are used to treat companion animals, ruminants, swine, and poultry [16]. Tulathromycin, the macrolide used in our study, was approved by the FDA in 2005 to treat and control respiratory disease in cattle and swine [16].
Antimicrobial treatment is important for treatment and prevention of infections by pathogenic bacteria, but may also select for resistant strains and increase the prevalence of resistance genes that may be transferred to other bacterial strains [1,18]. The majority of antimicrobial resistance studies have focused on the effect of antibiotic treatment on individual bacterial strains [19,20]. Next Generation Sequencing is a relatively novel technique that is helping expand our knowledge of the impact of antibiotic treatment on the gut microbiome and the prevalence of bacterial genes in diverse bacterial populations [5,21,22]. There is limited published information on the impacts of antimicrobial drug treatment on microbial genetic function in dairy calves, leaving a knowledge gap that limits the thorough evaluation of the impact of the use of wide spectrum antimicrobial drugs in young calves.
The objective of this study was to longitudinally characterize the effect of enrofloxacin or tulathromycin metaphylactic treatment on the fecal microbiome and microbial genetic function of preweaned dairy calves, focusing on functions related to antimicrobial resistance.

Ethics statement
The research protocol was reviewed and approved by the Institutional Animal Care and Use Committee of Cornell University. The drug administration to calves housed on the commercial dairy farm and the collection of fecal samples was authorized by the farm owner, who was aware of all experimental procedures. The herd veterinarian determined the need for metaphylactic treatment and the drugs were administered according to label directions.

Farm management
The study was conducted from May to November 2015 at a commercial dairy farm that milked 1,200 Holstein cows near Ithaca, New York, USA. The farm was selected because preventive antimicrobial treatment was indicated based on a prior history of calfhood respiratory disease identified by the herd veterinarian, as well as conditions that resulted in calves being at high risk of developing bovine respiratory disease (BRD), warranting use of antibiotic for control of BRD as per drug label. Routine calf management was performed by farm employees. Newborn calves received four quarts of colostrum within the first hours of life and were bottle fed pasteurized non-saleable milk twice a day until weaning at approximately 56 days of age. Water was offered ad libitum in buckets. Heifer calves were kept in individual hutches or pens with sawdust bedding during the preweaning period (first 8 weeks of life) and then moved to group pens. Health-related events (e.g. otitis, pneumonia and diarrhea) were treated as needed by farm employees and recorded in the farm's herd management software (Dairy-Comp 305; Valley Ag Software, Tulare, CA, USA). The research group obtained these health-related and treatment records from the herd management software and from the farm's drug use notebook to increase data accuracy.
A randomized field trial study design was used. Calves 2 to 3 weeks old, which had not been treated with antimicrobials before enrollment, were randomly allocated within block (enrollment day) to one of three study groups: control (CON, no antimicrobial), a single treatment of enrofloxacin at label dosages (ENR, Baytril 100, Bayer Corp. Agricultural Division, Shawnee Mission, KS, USA) or a single treatment tulathromycin at label dosages (TUL, Draxxin, Zoetis, Floham Park, NJ, USA). A total of 84 heifers were enrolled in the trial: 26 allocated to CON, 28 to ENR and 27 to TUL, distributed in 6 cohorts (6 enrollment days). Fecal samples were collected directly from the rectum of individual calves starting before the administration of the antimicrobial treatment on enrollment day (Day 0) and on 4, 14, 56, and 112 days after enrollment. Fecal samples were transported in a cooler with ice packs to the laboratory at the Cornell campus in Ithaca, NY, where they were aliquoted and stored at -20˚C until DNA extraction. For whole genome sequencing, a subset of 12 calves per study group was randomly selected after excluding calves with missing samples or incomplete data, and calves treated with antimicrobials or other drugs by farm personnel or the herd veterinarian after enrollment.

DNA extraction
Fecal samples were thawed and total DNA was extracted using the MoBio PowerSoil DNA isolation kit (MoBio Laboratories, Carlsbad, CA, USA) following the protocol previously used by Pereira et al. with few modifications [23]. Briefly, approximately 50 mg of feces was transferred to the PowerBead Tube, which was heat treated at 65˚C for 10 minutes and then 95˚C for 10 minutes. PowerBead Tubes were vortexed horizontally using the MoBio Vortex Adapter tube holder at a maximum speed for 10 minutes. The remaining DNA extraction procedure followed the standard protocol supplied by the kit manufacturer. The DNA concentration and purity were evaluated by optical density using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Rockland, DE, USA) at wavelengths of 230, 260 and 280 nm. The eluted DNA was stored at -20˚C until further processing.

Library preparation and illumina HiSeq sequencing
To increase measurement accuracy for the concentration for the DNA library, the starting DNA library was quantified using a fluorometric-based method, the Qubit dsDNA HS Assay system (Life Technologies Corporation, Carlsbad, CA, USA). An aliquot of each DNA sample was diluted to 0.2 ng/μl and used as input DNA to the Nextera XT DNA Library Preparation Kit (Illumina Inc. San Diego, CA), according to manufacturer's recommendations. As a quality control step after PCR-cleanup, a subset of libraries was run on an Agilent Technology Bioanalyzer to check the size distribution (bp). After library normalization and pooling, sequencing was performed at the Genomics Facility of the Cornell University Institute of Biotechnology with an Illumina HiSeq2500 platform, operating in Rapid Run Mode, paired-end 2 x 250 bp, one sample pool in both lanes of a two-lane flow cell.

MG-RAST analysis
The paired-end sample sequences were obtained from the Cornell University Institute of Biotechnology as two individual fastq files with one file per end. The two files were concatenated using Python and uploaded to the Metagenomics Rapid Annotation using Subsystem Technology (MG-RAST), which is an open-source server that analyzes sequencing data and provides taxonomic and genetic functional classification [24]. In summary, before submission for analyses at MG-RAST, standard screening and quality control procedures were selected, including removal of artificial replicate sequences, host (Bos Taurus) sequences and low quality sequences. The Refseq and SEED Subsystems databases were selected for taxonomic distributions (Phylum, Class, Order, Family, Genus and Species) and genetic functional assignments (4 levels of hierarchy), respectively, as outlined in Fig 1.

Statistical analyses
Phylogenetic and genetic functional relative abundances in percentages per sample were calculated using JMP Pro 11 (SAS Institute, Cary, NC). For our study, we selected linear discriminant analysis (LDA) because it allows the calculation of the best discriminating components based on foreknowledge of defined groups. As a clinical trial, this is of importance for our study given that groups are well defined, and our goals are focused on discriminating between these study groups. Phylogenetic taxa relative abundance variables were used as covariates in a linear discriminant analysis (LDA) models to evaluate changes in taxonomic composition (50 most abundant taxa) over time and by days after enrollment for each study group (P-value < 0.05). Stepwise selection was used and taxa with P-value < 0.1 were included in the LDA model and in the subsequent multivariable mixed logistic regression model, described below.
Data were displayed using the centroid (mean) and an ellipsis that represents a 95% confidence interval of the mean for each day.
A series of multivariable screening analyses using JMP Pro 11 was performed to determine which bacterial species and functions were most important to differentiate between study groups. False discovery rate (FDR) was used to correct for multiple comparisons in screening analysis. Given the large sample size, a strict cutoff (FDR-probability � 0.001) was used to minimize Type-1 statistical errors, as previously described [25]. For each category, a model was fitted for the most prevalent taxa or genetic functions after initial screening using FDR. The independent variables study group (CON, ENR or TUL), sample day (4, 14, 56, and 112) and interactions were included as fixed effects in all models. Because values at enrollment diverged among the three groups, Day 0 measurements were included as a covariate in the model. The effects block and individual animals nested within block were controlled in the models as random effects.

Shotgun sequencing
Whole genome sequencing of 179 samples produced 228,587,989 sequences totaling 57,397,337,159 basepairs (bp). The sample length was standardized to 251 bp when merging paired-end samples. One sample (calf in enrofloxacin group, Day 4) with a low number of sequences was excluded. Sequencing data by sample are listed in S1 Table and more information is available at MG-RAST in project number 20043 (https://www.mg-rast.org/linkin.cgi? project=mgp20043).

Effect of antibiotic treatment on taxonomic composition
Firmicutes, Bacteroidetes, Actinobacteria and Proteobacteria were the four major phyla and accounted for 96.4% of the microbiota on average. The top 20 most relatively abundant phyla are presented in Fig 2 by sample day for each study group. Relative abundance is a proportion of a taxa or gene, and does not represent absolute numbers.
Change in the microbiota composition was observed at the phylum level using LDA, as shown by the three-dimensional canonical plot (Fig 3). In this figure, each ellipse, which indicates the 95% confidence region to contain the true mean of the group of variables, represents a sampling day. The position of the ellipses shows how the microbiota changed over time, revealing a greater shift in microbial composition observed from Day 14 to Day 112. When analyzing the same dataset by sample day, there was a difference at the phylum level on Days 0, with calves treated with enrofloxacin being more distant from the others, and on Day 112, when control calves were separated from the treated calves (S1 Fig). Discriminant analysis was performed at the subsequent phylogenetic levels. At the genus level, the differences were less pronounced compared to phylum level (S2 Fig). When looking at the effect of antibiotic treatment on each taxon individually, including fixed and random effects, there was no significant effect of study group at the phylum level. As already described, because of the difference in baseline value between study groups on Day 0, it was included in the mixed model as a covariate.
The relative abundance of Desulfovibrionales, an order within the Deltaproteobacteria class, which belongs to phylum Proteobacteria, was significantly different among groups, being higher on Day 56 in calves treated with enrofloxacin (P-value 0.015). A similar pattern was observed for the family Desulfovibrionaceae (P-value 0.010) and genus Desulfovibrio (P-value 0.009, Fig 4).
In the discriminant analysis, 16 unique genera were selected for the model (P-value < 0.1). When evaluated using the multivariable mixed logistic regression model, only Bacteroides, Blautia, Coprococcus and Desulfovibrio were significantly affected by study group and day (Pvalues 0.035, 0.019, 0.027 and 0.009, respectively; Fig 4).

Effect of antibiotic treatment on microbiota genetic function
As seen with microbial taxonomy, the linear discriminant analysis revealed a change in microbial genetic function over time, which is shown by the three-dimensional canonical plot from Level 1 in Fig 5 (P-value < 0.05). As with microbial composition, the ellipses show a greater change in microbial genetic function relative abundance as the animals became older.
At Level 1, the interaction between study group and day had a significant effect on the category "Virulence, Disease and Defense"; however, the difference was mainly on Day 0. When Day 0 was included in the model as a covariate, this significance was lost. Also at Level 1, there was a significant difference in the relative abundance of "Clustering-based subsystems" (Pvalue < 0.05), as shown in Fig 6. Despite differences observed over the first weeks, at Day 112 functions appeared to be at similar relative abundances among all study groups.
At Level 2, calves in the enrofloxacin group had a higher relative abundance of genes with the "Transposable elements" function, which is part of Level 1 category "Phages, Prophages, Transposable elements, Plasmids", on Day 4 (P-value 0.034) and then reached levels similar to the other two groups, as shown in Fig 7. The distribution of "Resistance to antibiotics and toxic compounds" (RATC) is depicted by the boxplots in Fig 8. There was a large variation in "RATC" relative abundance on Day 0, especially for the enrofloxacin group. Although calves that received enrofloxacin started with a higher relative abundance of "RATC" on Day 0, we did not observe a rise in AMR genes after treatment with enrofloxacin. According to the mixed model results, "RATC" was not statistically different among the three study groups (P-value 0.095, S3 Fig).
At level 3, within "RATC", we found no significant difference in the relative abundance of "Resistance to Fluoroquinolones" among study groups (P-value 0.19). Relative abundance of the "RATC" function increased over time for all study groups, including for calves that did not receive antibiotic treatment, and it was higher on Day 56 for calves treated with tulathromycin (Fig 9). Also within "RATC", the relative abundance of "Multidrug efflux pump in Campylobacter jejuni (CmeABC operon)" was significantly different among groups (P-value 0.015), with calves treated with enrofloxacin having higher relative abundance on Day 4 and then The independent variables study group (CON = control, ENR = enrofloxacin, TUL = tulathromycin), sample day (4, 14, 56, and 112) and interactions were included as fixed effects in all models. Day 0 was included as a covariate in the model. The effects block and individual animals nested within block were controlled in the models as random effects. Asterisks indicate significant differences (P-value � 0.05) between day and study group. Error bars represent the standard error of the least square mean.

Discussion
Our results support that metaphylactic treatment of preweaned calves with enrofloxacin or tulathromycin did not cause major changes in fecal microbiota composition and gene Relative abundance of genetic functional categories within "RATC" by sample day represented by least square means. The independent variables study group (CON = control, ENR = enrofloxacin, TUL = tulathromycin), sample day (4, 14, 56, and 112) and interactions were included as fixed effects in all models. Day 0 was included as a covariate in the model. The effects block and individual animals nested within block were controlled in the models as random effects. Asterisks indicate significant differences (P-value � 0.05) between day and study group. Error bars represent the standard error of the least square mean.
https://doi.org/10.1371/journal.pone.0219635.g009 functions related to antimicrobial resistance. Despite the fact that calves were not treated with antibiotics before enrollment and were randomly assigned to one of the three study groups, these groups started the trial with different microbiota composition and function, including higher relative abundances of "Virulence, Disease and Defense" and "RATC" functions in calves assigned to enrofloxacin treatment. Nevertheless, after controlling for this difference in the analysis, we did not observe a significant increase in the relative abundance of genes with these functions after treatment with antimicrobial drugs.
When lactating cows are treated with antimicrobial or other drugs due to illnesses, their milk is withheld from sale because it may contain drug residues. This non-saleable milk, also called waste milk, is often fed to calves. Antibiotic residues found in waste milk from 34 farms in New York state were mainly β-lactams, tetracycline and sulfamethazine [26]. Enrofloxacin is approved only for non-lactating female dairy cattle less than 20 months of age, and extralabel use of this drug is strictly prohibited; therefore, there should not be residues from this drug class in milk fed to calves. Additionally, enrofloxacin was not used at the study farm as a treatment for calf pneumonia before the study began.
In our study, calves of all groups were fed pooled pasteurized non-saleable milk twice a day, and this could have affected the microbiota of calves. We continued to use this farm management because it is representative management practice on dairy farms. As an effect on study outcomes, the same milk was fed to different study groups, and because of that, we do no perceive it as a factor that would result in significant difference between study groups. We could however have the hypothesis the treatment of calves with enrofloxacin or tulathromycin in calves not fed waste milk could have resulted in different microbiota outcomes [25], however our study aimed to generate data of validity to common dairy practices.
Calves could have acquired drug resistant bacteria and resistant genes from their mothers during birth, from the waste milk used as feed [25,27] or from the environment [8,28], including cross-contamination through farm workers and feeding equipment [29]. Studies have found AMR genes in animals raised without antimicrobials as well [30,31]. Auffret et al. found differences in AMR genes depending on the diet (forage versus concentrate) fed to beef cattle not treated with antibiotics [31].
More than 20 years ago, Brown et al. discussed the increase in resistance to fluoroquinolones over time and emphasized the need to use these drugs judiciously in human and veterinary medicine [32]. Cummings et al. alerted to an increase in enrofloxacin resistance in bovine E.coli isolates in the northeastern United States from 2004 to 2011 [33]. In a controlled trial evaluating effects of low concentrations of drug residues in waste milk, "Resistance to fluoroquinolones" was found on the first day of life and increased over the following 6 weeks, even in calves that did not receive any parenteral antimicrobial or drug residue in their milk [25]. In our study, "Resistance to fluoroquinolones" was present even before the use of enrofloxacin on Day 0, and increased over time, regardless of study group.
Treatment with antimicrobial drugs did not have a significant effect on "Virulence, Disease and Defense" and "RATC", as we hypothesized. On the other hand, calves treated with enrofloxacin had higher relative abundance of "Transposable elements" on Day 4. Also known as transposons, these mobile genes can transfer functions within genomes and cause mutations [34]. The enrofloxacin group also had higher relative abundances of "Clustering-based subsystems" and lower relative abundances of "Multidrug efflux pump in Campylobacter jejuni (CmeABC operon)" on Day 14; however, it is not clear what these changes may represent.
Changes in taxonomic composition over time are expected as the animal grows, starts eating solids and the microbiome becomes more diverse. The age group of 2 to 3 weeks old was selected because of the higher risk for development of BRD at this stage of life, and enrofloxacin and tulathromycin are commonly used for treatment of heifers with BRD [11,35]. The parenteral use of these two antimicrobials did not cause significant deviations in the relative abundances of the four major phyla: Firmicutes, Bacteroidetes, Actinobacteria and Proteobacteria. These results are in agreement with other studies that evaluated fecal microbiome through shotgun and 16s rRNA sequencing [23,25,36,37], but in different proportions.
Bacteroides is a frequently studied in humans, Gram-negative, anaerobic bacterial genus, composed of several species and strains [38,39]. They are abundant in the gut, but can be opportunistically pathogenic such as during trauma and post-surgical infections, [38,39]. Bacteroides can develop AMR to many antimicrobials, as extensively reviewed by Wexler et al. [40]. An increase in resistance to fluoroquinolones has been reported among B. fragilis group strains from humans in the United States (moxifloxacin) [41] and in Spain (moxifloxacin and trovafloxacin) [42]. Despite having more Bacteroides at Day 0, calves treated with enrofloxacin did not show a major increase in the relative abundance of genes related to "RATC", while non-treated calves had an increase in Bacteroides relative abundance over time and in "RATC" as well.
Calves in the enrofloxacin group had an increase on Day 56 in the genus Desulfovibrio, its family and order. Desulfovibrio spp. are sulfate-reducing bacteria that produce hydrogen sulfide, which is both important for cell physiology and toxic to epithelial cells [43]. They are also opportunistic pathogens found in the environment [44] and GI tract of humans and other animals [45,46]. Because they are difficult to culture, there is a lack of information about this Gram-negative anaerobe, including antimicrobial susceptibility data. In a study with 23 Desulfovibrio isolates from four different species cultured from human specimens, D. fairfieldensis strains showed resistance to β-lactams and the fluoroquinolones ciprofloxacin and levofloxacin [47]. Culture-independent studies will bring more knowledge to these pathogens, but culturebased methods continue to be important for antimicrobial susceptibility surveillance.
The genus Blautia has been associated with human gut health, being reduced in patients with colorectal cancer [48]. Together with other carbohydrate-utilizing bacteria, higher relative abundances of Blautia and Coprococcus were found by Song et al. in the gut of 3-week old dairy calves, which coincides with rumen development, solid intake and greater concentration of short-chain fatty acids [49]. These two genera were temporarily increased in calves that received enrofloxacin in our study, which may be considered a good, unpredicted side effect.
Pitfalls of our study include the fact that, although calves were randomly enrolled in study groups, initial microbial composition and genetic function still differed between calves in different study groups at Day 0. The genetic functional analysis can only show the genes that were present in those samples, but actual physiological measurements were not made in this part of the study. Fecal samples were also plated and aerobically cultured for antimicrobial resistance profiles of E. coli isolates for another study. This study was performed on a single farm, which did not allow comparisons of effects across herds. Nevertheless, the calf management practices of the farm used in the study are typical for many other commercial dairy farms in New York State.
In our study, enrofloxacin or tulathromycin had minimal impacts on the microbial composition and genetic functional microbiota of calves over the study period (112 days) when used to prevent and control BRD. However, these results do not support the indiscriminate use of antibiotics. More studies are required to investigate the gut microbiome on different days post antimicrobial treatment and other scenarios. Finally, the authors would like to emphasize that antimicrobials should be used prudently in medicine and agriculture.
Supporting information S1 Table. Sequencing data by sample. More information is available at MG-RAST in project 20043 (https://www.mg-rast.org/linkin.cgi?project=mgp20043). (PDF) S1 Fig. Linear discriminant analysis of the 50 most abundant phyla. Changes for each study group on Day 0 and Day 112 after enrollment (P-value < 0.05). An ellipse indicates the 95% confidence region to contain the true mean of the variable (group). CON = control, ENR = enrofloxacin, TUL = tulathromycin. (TIF) S2 Fig. Linear discriminant analysis of the 50 most abundant genera. A) Changes in genera composition over time (P-value < 0.05). B) Changes for each study group by days after enrollment. An ellipse indicates the 95% confidence region to contain the true mean of the variable (day or group). CON = control, ENR = enrofloxacin, TUL = tulathromycin. (TIF) S3 Fig. Relative abundance of genes with "Resistance to antibiotics and toxic compounds" function by sample day represented by least square means. The independent variables study group (CON = control, ENR = enrofloxacin, TUL = tulathromycin), sample day (4, 14, 56, and 112) and interactions were included as fixed effects in all models. Day 0 was included as a covariate in the model. The effects block and individual animals nested within block were controlled in the models as random effects. Asterisks indicate significant differences (P-value � 0.05) between day and study group. Error bars represent the standard error of the least square mean. (TIF)