Skip to main content
  • Research article
  • Open access
  • Published:

Dose-dependent impact of oxytetracycline on the veal calf microbiome and resistome

Abstract

Background

Antibiotic therapy is commonly used in animal agriculture. Antibiotics excreted by the animals can contaminate farming environments, resulting in long term exposure of animals to sub-inhibitory levels of antibiotics. Little is known on the effect of this exposure on antibiotic resistance. In this study, we aimed to investigate the long term effects of sub-inhibitory levels of antibiotics on the gut microbiota composition and resistome of veal calves in vivo.

Forty-two veal calves were randomly assigned to three groups. The first group (OTC-high) received therapeutic oral dosages of 1ā€‰g oxytetracycline (OTC), twice per day, during 5 days. The second group (OTC-low) received an oral dose of OTC of 100ā€“200ā€‰Ī¼g per day during 7 weeks, mimicking animal exposure to environmental contamination. The third group (CTR) did not receive OTC, serving as unexposed control. Antibiotic residue levels were determined over time. The temporal effects on the gut microbiota and antibiotic resistance gene abundance was analysed by metagenomic sequencing.

Results

In the therapeutic group, OTC levels exceeded MIC values. The low group remained at sub-inhibitory levels. The control group did not reach any significant OTC levels. 16S rRNA gene-based analysis revealed significant changes in the calf gut microbiota. Time-related changes accounted for most of the variation in the sequence data. Therapeutic application of OTC had transient effect, significantly impacting gut microbiota composition between day 0 and day 2. By metagenomic sequence analysis we identified six antibiotic resistance genes representing three gene classes (tetM, floR and mel) that differed in relative abundance between any of the intervention groups and the control. qPCR was used to validate observations made by metagenomic sequencing, revealing a peak of tetM abundance at day 28ā€“35 in the OTC-high group. No increase in resistance genes abundance was seen in the OTC-low group.

Conclusions

Under the conditions tested, sub-therapeutic administration of OTC did not result in increased tetM resistance levels as observed in the therapeutic group.

Background

Antibiotic therapy is commonly used in animal agriculture both to treat ill animals and as a form of metaphylaxis, to prevent the development and spreading of infections in high-risk conditions. Low doses of antibiotics can also be added to animal feed to promote growth, although this practice has been banned in the European Union since 2006 [1]. The use of antibiotics can lead to the selection of resistant strains in livestock, shedding in the food chain and posing a risk for food safety and human health [2]. Furthermore, resistant strains, along with residual antibiotics, antibiotic metabolites and antibiotic resistance genes, eventually can reach waters and soils [3]. As a result, bacteria in these environments are frequently exposed to low concentrations of antibiotics, which can impact microbial ecosystems through community shifts and alterations of the resistome [4, 5].

Under laboratory conditions, it has been shown that selection for antibiotic resistance can occur even at very low antibiotic concentrations [6]. Selection mechanisms for resistance at sub-inhibitory concentrations are different from those active in the presence of lethal antibiotic concentrations [7]. Exposure to very low concentrations of antibiotics ā€“ even below the minimum inhibitory concentration (MIC) ā€“ elicits an adaptive resistance response [8], and can result in the presence of a wider range of resistant mutants, increased genotypic [9] and phenotypic [8] diversity, and higher rates of mechanisms that allow spreading of resistance, such as horizontal gene transfer (HGT) [10]. It has also been shown that resistance mutations in bacteria exposed to sub-MIC antibiotic concentrations are different from those of bacteria exposed to high concentrations [11], supporting the presence of alternative resistance mechanisms.

Many aspects of the development of antibiotic resistance in complex microbial ecosystems encountered in vivo in the presence of low concentrations of antibiotics still need to be elucidated. While pure culture laboratory studies allow detailed analysis of the impact of antibiotics under well controlled experimental conditions, such studies lack the complex ecological interactions encountered in vivo. Competition between microorganisms, host-microbe interactions and the impact of environmental conditions (including food/feed intake) exert ecological pressure, and may influence the net effect of antibiotics, especially at low concentrations.

The main objective of this study was to assess the effects of a long term exposure to a low dose of antibiotics on the veal calf gut microbiome in vivo. We performed an intervention study where we administered oxytetracycline to calves, either at a high therapeutic dose (oral therapy for 5 days), or at a low dose, to mimic exposure to environmental contamination. We then used a combination of methods (antibiotic analysis by the UPLC-MS/MS, 16S ribosomal sequencing, shotgun sequencing, quantitative PCR) to determine the effects of the interventions on the composition of the fecal microbiota of the calf and on the presence of antibiotic resistance genes therein.

Methods

Experimental design: Animals, treatment and sample collection

Farm-level studies were performed between 9ā€‰āˆ’ā€‰6-2015 and 28-7-2015 accordance with the Dutch Law on Animal Health and Welfare. Sixty healthy male Friesian Holstein calves (between 2 and 4ā€‰weeks of age), bodyweight approx. 50ā€‰kg, were collected from different dairy farms in in the regions Nordrhein-Westfalen, Rheinland-Pfalz, Saarland, ThĆ¼ringen, Sachsen and Hessen (Germany) and transported via a calf-collecting-center in Nordrhein-Westfalen (Germany) to a veal farm where they were housed individually in indoor pens in a clean, newly built barn. Animals were placed under regular agricultural care according to the regulations of the Integrated Chain Management System (IKB) and regular clinical veterinary practices. Following a 1 week period during which the calves adapted to the new environment and were fed according to routine, the intervention was initiated. Of the 60 calves, 42 calves that were shown to have fecal antibiotic levels upon arrival below 100ā€‰Ī¼g/kg feces were selected and assigned to three intervention groups. Each intervention group consisted of 14 animals, housed in such a way as to separate the control and low dose group from the therapeutic level intervention group. In addition, contact between calves was only possible with neighboring animals of the same experimental group (Additional file 1: Figure S2). Logistics in the barn were organized in such a manner that it would minimize carriage of antibiotic residues within the barn. In addition all materials required for animal treatment were antibiotic free and were available in separate color coded sets for each of the three groups. The first group (OTC-high) received an oral dose of 1ā€‰g OTC in two liters Calf Milk Replacer (CMR), administered twice a day for 5 days. The second group (OTC-low) received an oral dose of OTC of 25ā€‰Ī¼g/L CMR daily during 7 weeks (increasing from 100ā€‰Ī¼g/day at the start of the intervention until 200ā€‰Ī¼g/day at the end), mimicking animal exposure to environmental contamination [12, 13]The third group (control) did not receive OTC, serving as a unexposed control (Additional file 2: Figure S1). The three groups are referred to in the paper as OTC-high, OTC-low and CTR respectively. Calves were fed antibiotic-free CMR twice daily throughout the trial. The amount of CMR increased from four liters per day at day 0 to eight liters per day at day 42. The CMR was supplemented with etheric oil solution (Bronch-Arom F) and linseed oil. The calves also received roughage (starting with 50ā€‰g/d, rising to 450ā€‰g/d) and a homeopathic preparation throughout the trial, and they were on individual basis supplemented with Globatan, vitamin E and Vitamix.

On three occasions, individual treatment using non-antibiotic supplements did not sufficiently relieve respiratory symptoms. An antibiotic treatment was therefore required upon indication by the responsible veterinarian. In accordance with the study protocol, non-tetracycline antibiotics were given to all animals in order to avoid any bias between the study groups. Three days prior to the start of the study, all calves were given 600ā€‰mg tilmicosin (Milbosin, Dechra). At day 5, 1500ā€‰mg Florfenicol (NuFloR, MSD) was administered. At day 21, tilmicosin (Milbosin) was again given to all animals at a dose of 900ā€‰mg. Milbosin was given by subcutaneous injection. NuFlor treatment was done by intramuscular injection. Treatment was successful and relieved animals from respiratory symptoms and prevented further infection among the flock.

Fecal samples were collected within 24ā€‰h after arrival (day āˆ’ā€‰6), at the start of the intervention (day 0), and at days 2, 6, 14, 21, 28, 35 and 42. The fecal samples were cooled upon collection, and shipped to the laboratory where they were aliquoted and snap frozen in liquid nitrogen and stored at āˆ’ā€‰80ā€‰Ā°C. The time between sampling and freezing was at most 4ā€‰h.

Antibiotic residue analysis

After thawing, homogenized feces (1ā€‰Ā±ā€‰0.05ā€‰g) were weighed in a 50ā€‰mL polypropylene test tube. Samples were spiked with an internal standard of antibiotic isotope variants as appropriate and mixed. The antibiotics were extracted from the feces using 10ā€‰mL methanol:acetronitril:0.1ā€‰M EDTA-McIlvainbuffer (pHā€‰4.0), (30:20:50 v/v %) and 15ā€‰min shaking. After extraction the samples were centrifugated at 3800 x g for 15ā€‰min. The organic phase of the supernatant was evaporated under a stream of nitrogen at a temperature not exceeding 45ā€‰Ā°C. After this the extraction was repeated, and the two supernatants were combined, 27ā€‰mL 0.2ā€‰M Phosphate buffer (pHā€‰6.5) was added. Further clean-up was performed on Oasis HLB cartridges (200ā€‰mg, 6ā€‰cc, WATERS). 5ā€‰mL of the supernatant was passed through the columns. After elution the extracts were evaporated to dryness under a stream of nitrogen at a temperature not exceeding 45ā€‰Ā°C. The residue obtained was dissolved in 500ā€‰Ī¼L 0.1% HCOOH in H2O/MeOH (1:3 v/v %) for injection into the UPLC-MS/MS system. Analysis was performed on a Sciex 6500 QTRAP mass spectrometer (Sciex, Massachusetts, USA) connected to an Agilent 1290 Infinity (Agilent, California, USA) LC system. The mass spectrometer was operated in a positive and negative electrospray ionisation (ESI+ and ESIāˆ’) mode. The desolvation temperature was 400ā€‰Ā°C and the cone temperature was 120ā€‰Ā°C. Over 100 analytes were analysed in scheduled multiple reaction monitoring (SMRM) mode. Two transitions per analyte were created. The chromatography was performed on an Acquity UPLC BEH C18, 2.1x100mm, 1.7ā€‰Ī¼m column (WATERS, Massachusetts, USA). A gradient containing 0.1% formic acid (A) and 0.1% formic acid (HCOOH) in MeOH (B) was applied. The gradient went from 95% A stepwise to 100% B. The runtime for each injection was 12.5ā€‰min, the flow 0.4ā€‰mLā€‰mināˆ’ā€‰1 and the injection volume was 5ā€‰Ī¼L. Between runs, calibration samples were run composed of a mixture of 46 antibiotics.

DNA extraction and qPCR

For DNA isolation, fecal samples were thawed on ice and lysed by bead beating (Mini-BeadBeater-24, Biospec Products Bartlesville, USA) for 2ā€‰min at 2800 oscillations minuteāˆ’ā€‰1 in the presence of 300ā€‰Ī¼l of lysis buffer (Mag Mini DNA Isolation Kit, LGC ltd, UK), 500ā€‰Ī¼L zirconium beads (0.1ā€‰mm; BioSpec products, Bartlesville, OK, USA) and 500ā€‰Ī¼L phenol saturated with 10ā€‰mM Tris-HCl and 1ā€‰mM EDTA pHā€‰8.0 (Sigma). After centrifugation DNA was extracted using the Mag Mini DNA Isolation Kit (LGC ltd, UK) in accordance to the manufacturers recommendations. DNA quality was assessed by routine gel electrophoresis as well as by capillary electrophoresis on the Fragment Analyzer (Advanced Analytical, Heidelberg, Germany). Quantitative PCR (qPCR) primers and probes are listed in Additional file 3: Table S1. qPCR was performed using RT PCR master mix (Diagenode, Seraing, B) on an Applied Biosystems 7500 RT PCR system.

High-throughput sequencing

For 16S rDNA amplicon sequencing of the V4 hypervariable region, 1ā€‰ng of DNA was amplified as described by Kozich et al. [14] with the exception that 33ā€‰cycles were used instead of 35, using F515/R806 primers [15]. As control, mock samples comprising a mixture of 24 pure culture isolates, blanco extraction controls and two pooled fecal samples were included in each batch. Primers included the Illumina adapters and a unique 8-nt sample index sequence key [15]. The amount of DNA per sample was quantified using the Quant-iTā„¢ PicoGreenĀ® dsDNA Assay Kit (Thermo Fisher Scientific). The amplicon libraries were pooled in equimolar amounts and purified using the IllustraTM GFXTM PCR DNA and Gel Band Purification Kit (GE Healthcare, Eindhoven, The Netherlands). Amplicon quality and size was analysed on the Fragment Analyzer (Advanced Analytical). Paired-end sequencing of amplicons was conducted in five separate runs on the Illumina MiSeq platform (Illumina, Eindhoven, The Netherlands). The sequence data was processed with Mothur v.1.31.2 [16], in line with the mothur MiSeq SOP [14]. Sequences were grouped in operational taxonomic units by Minimal Entropy Decomposition (MED) using a minimum substantive abundance value (āˆ’M) of 500 [17]. Taxonomy was then assigned by querying the representative sequence of each oligotype against the SILVA database (release 132) [18].

Metagenomic shotgun sequencing (125ā€‰nt paired end sequencing) was performed at Baseclear (Leiden, The Netherlands) on the Hiseq2500. After quality filtering and removal of host genes, metagenomic sequence reads were mapped against the antibiotic resistance gene database CARD [19] (libraries AR, AT and ABS) using mapping software bowtie2 [20](applying methodā€‰=ā€‰āˆ’-very-fast-local). This allowed enumeration of reads/resistome element in each sample. Reads were allowed to map on multiple resistome genes (multiple mapping) in order to obtain maximal sensitivity to detect resistome elements. Mapping frequencies were normalized by using the total number of reads per fecal sample. An abundance filter was applied in which each resistance gene was required to have at least 3 sequence reads in at least four animals in at least one group at one timepoint. Metagenomic shotgun sequences were mapped against the non-redundant protein database through BLASTX using DIAMOND (V.0.9.9) [21]. Mapping results were parsed in MEGAN (v. 6) [22] and linked to KEGG functional modules (KEGG reference library december2017) [23].

Statistical analysis of 16S rRNA amplicon and metagenomic sequencing data

All statistical analyses on the 16S sequencing data were performed using R version 3.3.2 [24]. The overview of all sample groups in Fig.Ā 2 was creating using classical multidimensional scaling as implemented in R ā€˜baseā€™. A multidimensional scaling procedure was applied to a Bray-Curtis dissimilarity matrix created using the ā€˜distā€™ function of the ā€˜proxyā€™ package [25]. All distance displayed in the boxplots of Fig.Ā 3 were extracted from the same distance matrix as the one used for the multidimensional scaling procedure. Statistics on the distances shown in Fig. 2 were performed using linear mixed models with Tukey HSD post-hoc testing (no P-value correction), as implemented in the packages ā€˜lme4ā€™ and ā€˜multcompā€™ respectively [26, 27]. For the linear mixed model, the subject was used a random factor. Multivariate statistics on the 16S sequencing data was performed using PERMANOVA and constrained analysis of principal coordinates as implemented in the ā€˜veganā€™ package [28]. Shannon diversity was calculated using the ā€˜veganā€™ package. The diversity measures were transformed to deal with non-normality using Box-Cox transformation (as implemented in the ā€˜caretā€™ package [29]) and subsequently analyzed for statistical differences using linear mixed models with the subject as a random effect.

The KEGG module count data derived from the metagenomic shotgun data were analysed using DESeq2 [30], comparing differences in relative abundance of KEGG modules between groups at day 6 and day 42.

Statistical analysis of resistome

For each timepoint, the abundance of each gene as a function of the experimental group was modeled using a Poisson regression model with overdispersion [31]. The abundance was normalized for the total number of reads in each sample. The model specification is as follows:

$$ {Y}_{ij}\sim P\left({\mu}_{ij},\varPhi \right), $$

where Y denotes the abundance, i denotes the experimental group (1ā€‰=ā€‰control, 2ā€‰=ā€‰high, 3ā€‰=ā€‰low), j denotes the animal in a group (jā€‰=ā€‰1ā€¦9), is the true value of the abundance, P denotes a Poisson distribution with overdispersion, and Ī¦ denotes the dispersion parameter.

The true value of the abundance is related to the experimental conditions by

$$ \log {\mu}_{ij}=C+{\varDelta}_i+\log \raisebox{1ex}{${h}_{ij}$}\!\left/ \!\raisebox{-1ex}{$H$}\right. $$

where C is the level of control group, Ī”i is the increase of group i compared to the control, hij denotes the total number of reads of animal j in group i, and H denotes the largest number of reads in all 81 samples.

For each time point and gene, the Poisson regression model results in a P value that assesses whether the experimental groups differ in abundance. We applied a multiple test correction on all the P values [32]. The genes whose FDR-corrected P values for the test on group differences were smaller than 10% for at least one time point were studied further with t tests comparing the respective intervention groups with the control group.

Data deposition

The datasets supporting the conclusions of this article are available in the Sequence Read Archive (SRA) repository, project ID: PRJNA511412; https://www.ncbi.nlm.nih.gov/bioproject/PRJNA511412.

Results

Intervention

To investigate the effects of therapeutic and environmental oxytetracycline exposure on the calf fecal microbiome and resistome, an intervention study was performed. After a run-in period of 7 days, 42 calves of 2ā€“4ā€‰weeks of age were randomly assigned to three groups. The first group (OTC-high) received therapeutic oral dosages OTC, of two grams daily during 5 days. The second group (OTC-low) received a daily oral dose of 100ā€‰Ī¼g/day and increasing to 200ā€‰Ī¼g/day during 7 weeks, mimicking animal exposure levels to environmental contamination, as observed in previous studies [5, 12, 13, 33,34,35,36]. The third group (CTR) did not receive OTC, serving as a unexposed control. Animals were housed in individual pens, isolating the control and low dose group from the therapeutic level intervention group.

OTC levels discriminate between the intervention and control groups

Fecal antibiotic residue levels were determined by UPLC-MS/MS techniques. OTC, tetracycline, tilmicosin, tylosin and florfenicol were the only antibiotics detected to levels greater than 5ā€‰Ī¼g/kg in the fecal samples taken during the intervention period. The low levels and time of occurrence for tetracycline and tylosin were in correspondence with their likely origin as contaminants in the veterinary medical products used of OTC and tilmicosin administration. The levels of fecal OTC clearly discriminated the three groups (Fig.Ā 1). The therapeutic intervention group (OTC-high) showed a marked increase in fecal OTC levels immediately upon initiation of the treatment which declined to lower levels during the intervention. The median fecal OTC levels in the OTC-high group increased at day 2 to 335ā€‰mg/kg, and 682ā€‰mg/kg at day 6. At subsequent timepoints, the fecal OTC levels declined by approximately 90% each week. For the group receiving a continuous oral dose of OTC of 25ā€‰Ī¼g/L calf milk replacer (CMR), the median fecal levels of OTC were highest at day 6 at 111ā€‰Ī¼g/kg, declining to 8ā€‰Ī¼g/kg at day 42. Measurable fecal levels of OTC were detected in the control group, with highest median levels of 13ā€‰Ī¼g/kg at day 6. OTC fecal levels declined to 3ā€‰Ī¼g/kg at day 14 andā€‰<ā€‰2ā€‰Ī¼g/kg at day 28. In following time points no detectable levels of OTC was measured in feces. A marked consistency was found in the fecal OTC levels between calves within each group. In those samples in which higher levels of OTC was detected, we also found low levels of tetracycline. Tilmicosin was found in all three study groups, showing highest levels at day 0 and day 28, in line with application at day āˆ’ā€‰3 and day 21. Low levels of tylosin were detected in those samples in which highest levels of tilmicosin was found. Florfenicol was found in low but detectable levels (<ā€‰5ā€‰Ī¼g/kg) at day 14, probably as a residual trace of Florfenicol application at day 5.

Fig. 1
figure 1

OTC levels (Ī¼g/kg) measured at different time points in fecal samples from calves from the control, the low-dose, and the high-dose groups

Time and not antibiotic levels related most of the observed variation in microbiota composition

The fecal microbiota composition was analyzed for all calves at all time points by 16S ribosomal amplicon sequencing. In total 23,8 million raw sequence reads were generated, averaging 43.017 (Ā± 19.742) sequences per sample. After subsampling at 8.734 sequences, 606 operational taxonomic units were generated by minimum entropy decomposition, retaining a total of 3.4 million processed sequencesĀ (Additional file 4: Figure S3). To visualize the microbiota community development for the different groups over time we used multidimensional scaling (MDS) (Fig. 2A). The non-metric MDS plot showed that the largest variation (nMDS-1) in the microbiota composition was related to a time-dependent development. After day 14, the fecal microbiota community structure appeared to have stabilized, resulting in overlapping centroids of samples taken between day 14 and day 42. The second largest variation (nMDS-2), appeared to relate to shifts in microbiota community structure during the run-in period between day āˆ’ā€‰6 and day 0, linked to lowered Bacteroides levels and increased levels of Lactobacillus, Bifidobacterium and Blautia species. From the MDS plot, we did not see a clear separation between the control and the antibiotic intervention groups. To visualize possible differences related to the intervention groups, separated from temporal changes, we used Canonical Analysis of Principal coordinates (CAP), constraining for ā€˜timeā€™ (x-axis) and ā€˜groupā€™ (y-axis) (Fig. 2B). In the analysis we excluded samples from day āˆ’ā€‰6 to avoid a bias related to changes during the run-in period. This analysis showed that 11.4% of the variance could be explained by temporal changes between baseline (day 0) and day 42, revealing a similar temporal pattern for all groups. The variance linked to ā€˜groupā€™ differences was 1.2%, and was mainly related to a separation of the control group and the two intervention groups (high and low dose) at later time points. Microbiota changes related to time, as identified by examining the coefficients of the CAP model defining the time axis suggested a decrease in Alloprevotella, Bifidobacterium, Faecalibacterium, and Streptococcus species and increasing levels in Peptococcus, Prevotella and Bacteroidales species with increasing age (Additional file 5: Figure S4A). Microbiota changes related to groups related to higher le,vels of Ruminococcus, Coprobacillus and Lachnospiraceae species and lower levels of Prevotella, Faecalibacterium and Blautia species in the control group compared to the low and high dose intervention groups (Additional file 5: Figure S4B). Bacterial diversity (Shannon diversity) showed an increase over time, but showed no apparent difference between groups upon completion of antibiotic administration in the high group (day 6), nor upon completion of the study (day 42) (Additional file 6: Figure S5).

Fig. 2
figure 2

a Multidimensional scaling (MDS) plot of microbial community dissimilarities based on 16S rRNA gene sequencing data. The plot shows the relationships between the intervention and control groups. Samples are shown as translucent dots; the annotated opaque circles represent the centroids of each of the sample groups at the respective time points. b Canonical Analysis of Principal coordinates (CAP) of microbiota composition data. Analysis was constraint for time (x-axis) and group (y-axis)

We used permutational multivariate analysis of variance (PERMANOVA) to examine the significance of differences in microbiota community composition between intervention groups at each time point (TableĀ 1). With the exception of day 0 and day 14, the control and the OTC-high group were significantly different at each time point (pĀ <ā€‰0.05 at days 2 and 6; pĀ <ā€‰0.005 at days 21, 28, 35 and 42). The control and the low group were significantly different at day 21 through 42 (pā€‰<ā€‰0.05). The low and the high groups were significantly different at day 2 and 6 (pā€‰<ā€‰0.05) and at day 42 (pā€‰<ā€‰0.005).

Table 1 P-values derived from PERMANOVA pairwise statistical comparisons of microbial community composition between groups at each time point. Grey boxes indicate pĀ <ā€‰0.05; grey boxes with bold text indicate pĀ <ā€‰0.005

We also applied PERMANOVA to assess the significance of differences in microbiota community composition within groups between subsequent timepoints (TableĀ 2). The analysis showed significant changes within the control group during the first 2 weeks (PĀ <ā€‰0.05) but not in the succeeding time points. For the OTC-high as well as the OTC-low dose group, significant changes were detected at each subsequent time interval throughout the study (Pā€‰<ā€‰0.05). The microbiota differences supporting the PERMANOVA models were in full agreement with the coefficients defined in the Canonical Analysis of Principal coordinates (CAP) analysis related to temporal changes. We quantified the change in fecal microbiome composition for the individual calves at sequential time points by calculating the Bray-Curtis distances (Fig. 3). This confirmed the relative large community shifts between day āˆ’ā€‰6 and day 0 and the stabilization of the microbiota composition with increasing age, as noted by the smaller distances between subsequent time points and smaller variation between calves. We found a significantly (Pā€‰<ā€‰0.05) larger Bray Curtis distance between day 0 and day 2 for calves in the OTC-high group compared to the other groups, indicating that within this timeframe a significantly larger change in the microbiome occurred in the OTC-high group compared to other groups. Between day 2 and day 6, the Bray Curtis dissimilarities of the OTC-high group appeared to be larger than in the other groups, but this was not statistically significant. No significant differences were detected in Bray Curtis distances between groups at subsequent time points.

Table 2 P-values derived from PERMANOVA pairwise statistical comparisons of microbial community composition within groups between subsequent time point. Grey boxes indicate pā€‰<ā€‰0.05; grey boxes with bold text indicate pā€‰<ā€‰0.005
Fig. 3
figure 3

Boxplots showing the Bray-Curtis dissimilarities between individual calves within each group, for pairs of subsequent time points. The boxes in the plot represent the interquartile ranges, the horizontal lines give the position of the medians, the vertical bars indicate the range. The dots indicate outliers

Functional differences discriminate both intervention groups from the control group

We performed metagenomic shotgun sequencing of feces samples taken at baseline (day 0), day 6 and day 42, to explore potential functional differences in the fecal microbiota between the intervention groups. To facilitate the comparison, we aggregated the shotgun reads to the level of functional gene orthologous, using the KEGG Orthology, and we used DESeq2 to reveal significant differences between intervention groups. Table S2 lists the orthologous groups with significantly different abundances between groups (pĀ <ā€‰0.01). Most functional differences were detected between groups at day 42. At day 6, one gene orthologue (K19172) was found to be enriched in the high group compared to the control and low group. This orthologue group was found to be depleted in the high group compared to the control and low groups at day 42. K19172 consists of DndE, DNA sulfur modification proteins [37]. No other functional differences were detected between groups at day 6. The most profound difference detected at day 42 was the higher abundance of orthologous gene groups K18231 and K11959 in the high group compared to both the control and the low group. K18231 showing the largest fold change difference comprises membrane transporters involved in resistance to macrolides. K11959 also comprises transporter genes, linked to the urea. We detected overlap in the functional differences between either intervention group (high or low) and the control. Of the gene orthologues that were significantly different between the high or the low group and the control, four (K00689, K03620, K03605 and K02404) were significantly different in both groups compared to the control. K00689, K03620 and K03605 showed lower abundance in both the low and high group compared to the control, while K02404 was enriched in both groups. K00689 consists of orthologues of dextransucrase, a glycosyltransferase involved in carbohydrate metabolism. Orthologue group K03620 contains a hydrogenase and cytochrome subunit. K03605 is a protease that, by cleaving off a 15 amino-acid peptide, is involved in the terminal processing of a metal-binding hydrogenase. K02404 is a protein necessary for the biosynthesis of the flagellum, and may be involved in translocation of the flagellum.

Antibiotic resistance genes are significantly more abundant in the high intervention group at day 42

Antibiotic resistance (AbR) genes in the calf fecal microbiota were identified by aligning the metagenomic sequences of the calf fecal samples taken at day 0, 6 and 42 against the CARD database, a curated database of antibiotic resistance genes. Of the 532 antibiotic resistance genes identified, seven genes showed a significant difference in their relative abundance between either intervention group and the control group in at least one time point. Of these seven genes, two AbR genes had lower and five AbR genes had higher abundance levels in one of the intervention groups compared to the control (Table S3). Within the group of five genes with increased abundance, two genes conferred resistance to florfenicol (flo/floR), two genes represented macrolide resistance genes (mel), and one conferred resistance to tetracycline (tetM). The seventh gene identified by statistical comparison was RNA polymerase beta-prime chain rpoB, part of the CARD database as target for rifamycin. The normalized read abundances of the mel, tetM and flo genes in the three groups at the three time points is plotted in Fig.Ā 4.

Fig. 4
figure 4

Boxplots showing normalized read abundance of the antibiotic resistance genes tetM flo and mel for each group at different time points. Gene abundance is based on metagenomic data. Antibiotic resistance genes were identified for having a significantly higher abundance in the metagenomic database in the treatment group compared to the control group

Compared to the control group, the tetM tetracycline normalized gene abundance was found to be significantly higher in the therapeutic OTC-high group (PĀ <ā€‰0.05). The median normalized gene counts at day 42 were for the control group 10.7, for the OTC-low group 5.3 while for the OTC-high group this was 32.9. It was noted that at day 0, generally the normalized tetM gene counts were higher compared to other timepoints analysed. Median normalized tetM gene count was 48 at baseline and no significant difference in abundance was detected between groups. At day 6, the median normalized tetM gene count was 11,5.

For the mel gene encoding a macrolide efflux pump, significantly (Pā€‰<ā€‰0.05) higher gene counts were found for the OTC-high group. This was most pronounced at day 42. While the median gene counts at baseline were low (between 0 and 2.6), at day 6 the median gene count for the OTC-high group was 51. At day 42, the median gene count was 10.3 for the control group, 5.3 for the OTC-low group, and 90 for the OTC-high group. Considerable variation in the normalized gene counts was observed between calves within the OTC-high group at day 42.

Flo and floR abundance in general was low. Significant higher levels were found at day 42 in the OTC-high and OTC-low group in compared to the control group. The median levels at day 42 were 0 for the control and 2.9 and 1.0 for the OTC-high and OTC-low group respectively.

qPCR analysis confirms differences in antibiotic resistance genes abundance between groups

To confirm the findings obtained in the resistome metagenomic analysis and to further analyse the dynamic changes in the abundances of the resistance genes identified, quantitative PCR was performed. For this, we retrieved the metagenomic reads mapping to the mel, tetM and floR resistance genes and aligned them with corresponding genes from the CARD database. Next, we designed dedicated primers and probes and performed qPCR analysis in all 432 samples collected throughout the intervention. The mel, tetM and floR abundancies were normalized for the total bacterial load as established by broad range 16S ribosomal PCR. The results of the qPCR are showed in Fig.Ā 5. qPCR results confirmed the elevated tetM and mel abundance in the group receiving therapeutic OTC levels compared to the control in sample acquired at day 42. The relative abundance of tetM appeared to be highest between day 28 and day 35 for the OTC-high group. The tetM levels did not increase in the samples taken from the OTC-low group. We did detect elevated levels in the samples taken at the start of the intervention, declining to relatively low levels at day 14. The macrolide efflux gene, mel, was found to have higher abundance in the therapeutic OTC-high group, compared to the control group. The OTC-low group did not show elevated mel levels. floR had extremely low levels, near to the qPCR detection limit. It showed a very modest increase at day 14 and remained elevated up until day 28.

Fig. 5
figure 5

Relative abundance levels of antibiotic resistance genes TetM (upper row), Mel (middle) and floR (lower row) determined by real-time qPCR for each group and time point. The boxes in the plot represent the interquartile ranges, the horizontal lines give the position of the medians, the whiskers indicate the minimum and the maximum values. The dots indicate outliers

Discussion

In this study we aimed to investigate the effects of long term exposure to a low dose of oxytetracycline (OTC) on the calf fecal microbiota and resistome in vivo in comparison to short term therapeutic OTC treatment and a non-treated control. During the study, calves either received an oral therapeutic dosage of OTC for 5ā€‰days, or a low oral dosage of OTC for 7 weeks. Fecal samples were collected at weekly intervals and analysed in a combinatory approach using 16S ribosomal sequencing, shotgun metagenomics, quantitative PCR and UPLC-MS/MS analysis. In the therapeutic group (OTC-high), the fecal OTC levels were at least tenfold higher than the levels established for OTC resistant isolates of 16ā€‰Ī¼g/ml as reported by the Clinical and Laboratory Standards Institute. In the OTC-low group, OTC levels were hundredfold below these values and were in line with levels detected in contaminated farm environments [12, 13, 33, 34, 36]. Measurable levels of OTC were detected temporarily in the control group despite measures to prevent this. These levels were at least ten-fold lower than those measured in the OTC-low group. The set up thus allowed us to investigate the effects of OTC at therapeutic and sub-therapeutic concentrations and to compare these effects to a group that experienced negligible levels of OTC. However, it should be noted that herd-level application of therapeutic levels of tilmicosin and florfenicol was required during the study for the treatment of respiratory infections.

We observed profound changes in the fecal microbial community structure during the run-in period, between day āˆ’ā€‰6 and day 0 (baseline), coinciding with a decrease of Bacteroides and an increase of Blautia, Bifidobacteria and Lactobacillus species (Fig. 2A). A change in diet may be the main reason for the observed drastic changes in community structure between days āˆ’ā€‰6 and 0. A more gradual change in fecal microbiota composition was observed between day 0 and day 21, during which Prevotella became the dominant genus. After day 21, the fecal microbiota community structure stabilized. Contrary to our expectations, the provision of therapeutic levels of OTC did not result in apparent shifts in calf gut community structure, nor in species diversity, when compared to the control and low-dose conditions. A significant, transient effect of the therapeutic oral dosage of OTC between day 0 and day 2 was determined, as indicated by a larger Bray-Curtis dissimilarity between sample of the same calf taken at these time points (Fig. 3). Whether oral dosage of OTC also impacts other parts of the gastro intestinal tract remains to be determined. Several studies have described the temporal changes in the calf fecal microbiome over time, and suggested ties with changes in the calf diet as well as physiological developments of the gastro intestinal tract [38, 39]. Meale et al. [40] indeed showed that a delayed transition from a high milk diet to an exclusively solid feed diet resulted in a more gradual shift in gut microbiota community structure and functionality. Few studies have investigated the effects of therapeutic application of OTC on the calf fecal microbial community. Oultram et al. [41] presented findings of OTC treatment in seven calves, as part of a prospective cohort study in calves. Although the group size of the study of Oultram et al. [41] was small, discriminant analysis of the five most abundant genera revealed differences in the prevalence of the genus Lactobacillus and increased species richness in the treated calves compared to the controls. Antibiotic treatment also related to differences in health status between the animals.

To examine the effects of antibiotic administration on selection for antibiotic resistance genes, we performed metagenomic shotgun sequencing of samples taken at baseline (day 0), at day 6 and at day 42. After mapping metagenomic reads to known antibiotic resistance genes, statistical analysis was used to identify significant differences in gene abundance in any of the intervention groups compared to the control (Fig. 4). The tetracycline resistance gene tetM, encoding a ribosome protection protein (RPP) was shown to occur in significantly higher abundance levels in therapeutic group (OTC-high) at day 42. This was not observed for the group receiving the low oral dosage of OTC. Quantitative PCR confirmed elevated levels of tetM in the therapeutic group and showed that highest levels were measured at days 28 and 35 (Fig. 5). qPCR revealed no selective enrichment of tetM in the group receiving the low OTC dosage. In the group receiving the oral therapeutic dosage of OTC, significantly higher levels of the macrolide resistance gene mel were observed. In the group receiving the low dosage of OTC, we did not observe an increase in the levels of mel. Quantitative PCR confirmed the observations and showed that elevated levels of mel in the oral therapeutic group occurred at all time points after day 6. This coincides with the higher tilmicosin levels throughout the study. Since tilmicosin was given to all calves in all three groups, the significantly higher levels of mel resistance in the therapeutic OTC groups may suggest added effects of tetracycline and tilmicosine or multidrug resistance. In various studies on Gram positive clinical isolates, associations have been reported on tetracycline and macrolide resistance [42, 43]. In Streptococcus pneumoniae, resistance to macrolides and tetracyclines has linked to a conjugative transposon carrying multiple resistance genes [44, 45]. Phage-mediated transduction of mefA and tetO has been demonstrated between Streptococcus species [46]. Multidrug resistant Clostridium perfringens isolates from water, soil and sewage were shown to be capable of the transfer of tetM, ermB, ermQ and the mel orthologue mefA genes to the Enterococcus faecalis recipient through independent conjugative processes [47]. Resistome analysis also revealed an significant increase in the abundance of floR but its levels remained relatively low.

It has been well established that therapeutic levels of antibiotics can lead to selection of resistant bacteria. The effects of low concentrations of antibiotics, as can be encountered in farming environments, are not well understood. In vitro studies using high sensitive competition assays revealed that selection of resistant bacteria can occur at concentrations of more than hundred-fold below the minimal inhibitory concentration (MIC). For example, competition growth studies with wild type and resistant Salmonella typhimurium revealed that the lowest concentration at which novo resistant mutants were selected (minimum selective concentration) was 15ā€‰Ī¼g/L for tetracycline [6]. Mathematical modelling next showed that competitive advantage would allow such resistant mutants to become dominant in the population of susceptible organisms. Exposure to sub-MIC antibiotic concentrations may therefore be an important mechanisms for the selection and spread of resistance [48]. However, results from studies on the effects of subtherapeutic levels of antibiotics in vivo have not been conclusive, likely reflecting differences in set-up and type of antibiotics used. The effects of low, subtherapeutic concentrations of antibiotics on the fecal microbiota may be more complex as multiple factors drive selection in this ecosystem, including diet, immunity, and physiological conditions. This hampers the direct translation of effects detected in vitro to conditions in vivo. In a study, preweaned calves were fed for 6 weeks with milk spiked with low concentration of an antibiotic mixture of ceftiofur (100ā€‰Ī¼g/L), penicillin G (5ā€‰Ī¼g/L), ampicillin (10ā€‰Ī¼g/L) and OTC (300ā€‰Ī¼g/L), resulting in higher proportions of antibiotic resistant Escherichia coli and a number of genus-level changes in the fecal microbiota composition [49, 50]. Another study, instead, examined the effects of therapeutic (1ā€‰g/calf/day for 14ā€‰days) and subtherapeutic (10ā€‰mg/calf/day for 42ā€‰days) levels of neomycin and OTC, and found no changes in the absolute (copies per gram of wet feces), nor the relative (copies per copy of the ribosomal 16S gene) abundance of antibiotic resistance genes corresponding to the tetracycline (tetC, tetG, tetW, and tetX), macrolide (ermB, ermF), and sulfonamide (sul1, sul2) classes of antibiotics, nor for the class I integron gene, intI1. In that study, the relative abundance of tetO was shown to be higher in the therapeutic group [51].

Conclusions

The aim of this study was to investigate the impact of a long term exposure to a low dose of oxytetracycline on the gut microbiome and resistome of calves in comparison to therapeutic treatment and a non-treated control. Based on 16S rDNA sequence analysis, we observed that antibiotic administration did not cause drastic changes in the composition of the gut microbiome of calves. Within the first 2ā€‰days from the start of the intervention, a significant difference was observed between the high dose group and the other groups (control and low dose of antibiotics). However, after day 2, no significant differences were observed between changes occurring in the different groups. Multidimensional scaling analysis showed that time rather than group explained the largest variation in the microbiota between day 0 and day 21, and that after day 21 the microbial community structure appears to have stabilized. Resistome analyses revealed elevated levels of the tetracycline resistance gene tetM in the group receiving the oral therapeutic dosage of OTC as well as increased levels of the macrolide efflux gene mel in both therapeutic groups. This was not observed in the group receiving the low oral dosage of OTC throughout the 7ā€‰weeks study. As in complex microbial communities, multiple factors can act simultaneously to apply selective forces on members of this ecosystem, the net result of subtherapeutic levels of antibiotics are difficult establish. In vitro studies are likely to fall short as they lack some of the key drivers in ecological interactions and probably overestimate the impact of selective agents. In this study, metagenomic sequencing allowed us to establish the dose-dependent impact of antibiotics on the fecal microbial community structure and antibiotic resistance gene abundance in vivo. The unbiased nature of the methodology may be valuable in our efforts to establish minimum selective concentration values for antibiotics.

Abbreviations

CMR:

Calf Milk Replacer

CTR:

Abbreviation refers to the group of animals receiving not receiving any antibiotics and serving as unexposed control

HGT:

Horizontal gene transfer

MDS:

Multidimensional scaling

MED:

Minimal Entropy Decomposition

MIC:

Minimum inhibitory concentration

OTC ā€“high:

Abbreviation refers to the group of animals receiving the oral therapeutic dose of oxytetracycline

OTC-low:

Abbreviation refers to the group of animals receiving the oral low dose of oxytetracycline

qPCR:

Quantitative polymerase chain reaction

UPLC-MS/MS:

Ultra performance liquid chromatography - tandem mass spectrometer

References

  1. 1831/2003/EC otEPatCotEU: Regulation (EC) No 1831/2003 of the European Parliament and of the Council of 22 September 2003. Off J Eur Union 2003, 4:L 268/229ā€“242.

  2. Alban L, Ellis-Iversen J, Andreasen M, Dahl J, SĆønksen UW. Assessment of the risk to public health due to use of antimicrobials in pigs ā€“ an example of pleuromutilins in Denmark. Front Veterinary Sci. 2017;4:74.

    ArticleĀ  Google ScholarĀ 

  3. Yeom JR, Yoon SU, Kim CG. Quantification of residual antibiotics in cow manure being spread over agricultural land and assessment of their behavioral effects on antibiotic resistant bacteria. Chemosphere. 2017;182:771ā€“80.

    ArticleĀ  CASĀ  Google ScholarĀ 

  4. Jia S, Zhang X-X, Miao Y, Zhao Y, Ye L, Li B, Zhang T. Fate of antibiotic resistance genes and their associations with bacterial community in livestock breeding wastewater and its receiving river water. Water Res. 2017;124:259ā€“68.

    ArticleĀ  CASĀ  Google ScholarĀ 

  5. He LY, Ying GG, Liu YS, Su HC, Chen J, Liu SS, Zhao JL. Discharge of swine wastes risks water quality and food safety: antibiotics and antibiotic resistance genes from swine sources to the receiving environments. Environ Int. 2016;92-93:210ā€“9.

    ArticleĀ  CASĀ  Google ScholarĀ 

  6. Gullberg E, Cao S, Berg OG, IlbƤck C, Sandegren L, Hughes D, Andersson DI. Selection of resistant bacteria at very low antibiotic concentrations. PLoS Pathog. 2011;7:1ā€“9.

    ArticleĀ  Google ScholarĀ 

  7. Andersson DI, Hughes D. Microbiological effects of sublethal levels of antibiotics. Nat Rev Microbiol. 2014;12:465ā€“78.

    ArticleĀ  CASĀ  Google ScholarĀ 

  8. Sandoval-Motta S, Aldana M. Adaptive resistance to antibiotics in bacteria: a systems biology perspective. Wiley Interdiscip Rev Syst Biol Med. 2016;8:253ā€“67.

    ArticleĀ  Google ScholarĀ 

  9. Mira PM, Meza JC, Nandipati A, Barlow M. Adaptive landscapes of resistance genes change as antibiotic concentrations change. Mol Biol Evol. 2015;32:2707ā€“15.

    ArticleĀ  CASĀ  Google ScholarĀ 

  10. Jutkina J, Rutgersson C, Flach CF, Joakim Larsson DG. An assay for determining minimal concentrations of antibiotics that drive horizontal transfer of resistance. Sci Total Environ. 2016;548-549:131ā€“8.

    ArticleĀ  CASĀ  Google ScholarĀ 

  11. Westhoff S, van Leeuwe TM, Qachach O, Zhang Z, van Wezel GP, Rozen DE. The evolution of no-cost resistance at sub-MIC concentrations of streptomycin in Streptomyces coelicolor. ISME J. 2017;11:1168ā€“78.

    ArticleĀ  CASĀ  Google ScholarĀ 

  12. Pikkemaat MG, Yassin H, H.J. van der Fels-Klerx JJ, Berendsen BJA: Antibiotic Residues and Resistance in the Environment. In., vol. 2016.009: RIKILT Wageningen UR; 2016: 32.

  13. Schijndel JV, Oosterwegel J, Liefers R, Schmitt H, Schilt R, Lahr J. Antibiotica in de bodem. Een pilotstudie. Gouda: Stiching Kennisontwikkeling en kennisoverdracht Bodem (SKB); 2009.

    Google ScholarĀ 

  14. Kozich JJ, Westcott SL, Baxter NT, Highlander SK, Schloss PD. Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform. Appl Environ Microbiol. 2013;79(17):5112ā€“20.

    ArticleĀ  CASĀ  Google ScholarĀ 

  15. Caporaso JG, Lauber CL, Walters WA, Berg-Lyons D, Lozupone CA, Turnbaugh PJ, Fierer N, Knight R. Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample. Proc Natl Acad Sci U S A. 2011;108(Suppl 1):4516ā€“22.

    ArticleĀ  CASĀ  Google ScholarĀ 

  16. Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, Lesniewski RA, Oakley BB, Parks DH, Robinson CJ, et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol. 2009;75:7537ā€“41.

    ArticleĀ  CASĀ  Google ScholarĀ 

  17. Eren AM, Morrison HG, Lescault PJ, Reveillaud J, Vineis JH, Sogin ML. Minimum entropy decomposition: unsupervised oligotyping for sensitive partitioning of high-throughput marker gene sequences. ISME J. 2015;9(4):968ā€“79.

    ArticleĀ  CASĀ  Google ScholarĀ 

  18. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, Peplies J, Glƶckner FO. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41:590ā€“6.

    ArticleĀ  Google ScholarĀ 

  19. McArthur AG, Waglechner N, Nizam F, Yan A, Azad MA, Baylay AJ, Bhullar K, Canova MJ, De Pascale G, Ejim L, et al. The comprehensive antibiotic resistance database. Antimicrob Agents Chemother. 2013;57:3348ā€“57.

    ArticleĀ  CASĀ  Google ScholarĀ 

  20. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9:357ā€“9.

    ArticleĀ  CASĀ  Google ScholarĀ 

  21. Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2014;12:59ā€“60.

    ArticleĀ  Google ScholarĀ 

  22. Huson D, Mitra S, Ruscheweyh H. Integrative analysis of environmental sequences using MEGAN4. Genome Res. 2011;21:1552ā€“60.

    ArticleĀ  CASĀ  Google ScholarĀ 

  23. Kanehisa M, Goto S, Sato Y, Kawashima M, Furumichi M, Tanabe M. Data, information, knowledge and principle: Back to metabolism in KEGG. Nucleic Acids Res. 2014;42:199ā€“205.

    ArticleĀ  Google ScholarĀ 

  24. R Core Team: R: a language and environment for statistical computing. In. R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org/.

  25. Meyer D, Buchta C. Proxy: distance and similarity measures. R package version 0.4-17. 2017. https://CRAN.R-project.org/package=proxy.

  26. Bates D, MƤchler M, Bolker B, Walker S. Fitting linear mixed-effects models using {lme4}. J Stat Softw. 2015;67:1ā€“48.

    ArticleĀ  Google ScholarĀ 

  27. Hothorn T, Bretz F, Westfall P. Simultaneous inference in general parametric models. Biom J. 2008;50:346ā€“63.

    ArticleĀ  Google ScholarĀ 

  28. Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, Minchin PR, O'Hara RB, Simpson GL, Solymos P et al: vegan: Community Ecology Package. 2017.

    Google ScholarĀ 

  29. Kuhn M, Wing J, Weston , Williams A, Keefer C, Engelhardt A, Cooper T, Mayer Z, R Core Team, Benesty M, Lescarbeau R, Ziem A, Scrucca L, Tang Y, Candan C, Hunt T. caret: Classification and Regression Training. R package version 6.0-71. 2017. https://CRAN.R-project.org/package=caret.

  30. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:1ā€“21.

    ArticleĀ  Google ScholarĀ 

  31. McCullagh P, Nelder JA. Generalized Linear Models, Second Edition. 2 edition. Boca Raton: Chapman and Hall/CRC; 1989. p.532.

  32. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. ā€‹Journal of the Royal Statistical Society. Series B (Methodological). 1995;57(1):289ā€“300.

  33. Hamscher G, Pawelzick HT, Sczesny S, Nau H, Hartung J. Antibiotics in dust originating from a pig-fattening farm: a new source of health hazard for farmers? Environ Health Perspect. 2003;111:1590ā€“4.

    ArticleĀ  CASĀ  Google ScholarĀ 

  34. Kim KR, Owens G, Kwon SI, So KH, Lee DB, Ok YS. Occurrence and environmental fate of veterinary antibiotics in the terrestrial environment. Water Air Soil Pollut. 2011;214:163ā€“74.

    ArticleĀ  CASĀ  Google ScholarĀ 

  35. MuƱoz M, Autenrieth R, Mendoza-sanchez I. Environmental occurrence of Oxytetracycline and the potential selection of antibiotic resistance in bacteria. In: XVI world water congress: 2017; Cancun. Mexico: IWRA; 2017. p. 1ā€“15.

    Google ScholarĀ 

  36. De Liguoro M, Cibin V, Capolongo F, Halling-Sorensen B, Montesissa C. Use of oxytetracycline and tylosin in intensive calf farming: evaluation of transfer to manure and soil. Chemosphere. 2003;52(1):203ā€“12.

    ArticleĀ  Google ScholarĀ 

  37. Hu W, Wang C, Liang J, Zhang T, Hu Z, Wang Z, Lan W, Li F, Wu H, Ding J, et al. Structural insights into DndE from <i>Escherichia coli</i> B7A involved in DNA phosphorothioation modification. Cell Res. 2012;22:1203ā€“6.

    ArticleĀ  CASĀ  Google ScholarĀ 

  38. Uyeno Y, Sekiguchi Y, Kamagata Y. rRNA-based analysis to monitor succession of faecal bacterial communities in Holstein calves. Lett Appl Microbiol. 2010;51:570ā€“7.

    ArticleĀ  CASĀ  Google ScholarĀ 

  39. Oikonomou G, Teixeira AGV, Foditsch C, Bicalho ML, Machado VS, Bicalho RC. Fecal Microbial Diversity in Pre-Weaned Dairy Calves as Described by Pyrosequencing of Metagenomic 16S rDNA. Associations of Faecalibacterium Species with Health and Growth. PLoS One. 2013;8(4):e63157.

  40. Meale SJ, Li SC, Azevedo P, Derakhshani H, DeVries TJ, Plaizier JC, Steele MA, Khafipour E. Weaning age influences the severity of gastrointestinal microbiome shifts in dairy calves. Sci Rep. 2017;7:198.

    ArticleĀ  CASĀ  Google ScholarĀ 

  41. Oultram J, Phipps E, Teixeira AGV, Foditsch C, Bicalho ML, Machado VS, Bicalho RC, Oikonomou G. Short communication: effects of antibiotics (oxytetracycline, florfenicol or tulathromycin) on neonatal calves' faecal microbial diversity. Veterinary Record. 2015;177:598.

    CASĀ  PubMedĀ  Google ScholarĀ 

  42. Culebras E, Rodriguez-Avial I, Betriu C, Redondo M, Picazo JJ. Macrolide and tetracycline resistance and molecular relationships of clinical strains of Streptococcus agalactiae. Antimicrob Agents Chemother. 2002;46(5):1574ā€“6.

    ArticleĀ  CASĀ  Google ScholarĀ 

  43. Rubio-Lopez V, Valdezate S, Alvarez D, Villalon P, Medina MJ, Salcedo C, Saez-Nieto JA. Molecular epidemiology, antimicrobial susceptibilities and resistance mechanisms of Streptococcus pyogenes isolates resistant to erythromycin and tetracycline in Spain (1994-2006). BMC Microbiol. 2012;12:215.

    ArticleĀ  Google ScholarĀ 

  44. Courvalin P, Carlier C. Transposable multiple antibiotic resistance in Streptococcus pneumoniae. Mol Gen Genet. 1986;205(2):291ā€“7.

    ArticleĀ  CASĀ  Google ScholarĀ 

  45. Shen X, Yang H, Yu S, Yao K, Wang Y, Yuan L, Yang Y. Macrolide-resistance mechanisms in Streptococcus pneumoniae isolates from Chinese children in association with genes of tetM and integrase of conjugative transposons 1545. Microb Drug Resist. 2008;14(2):155ā€“61.

    ArticleĀ  CASĀ  Google ScholarĀ 

  46. Giovanetti E, Brenciani A, Morroni G, Tiberi E, Pasquaroli S, Mingoia M, Varaldo PE. Transduction of the Streptococcus pyogenes bacteriophage Phim46.1, carrying resistance genes mef(A) and tet(O), to other Streptococcus species. Front Microbiol. 2014;5:746.

    PubMedĀ  Google ScholarĀ 

  47. Soge OO, Tivoli LD, Meschke JS, Roberts MC. A conjugative macrolide resistance gene, mef(a), in environmental Clostridium perfringens carrying multiple macrolide and/or tetracycline resistance genes. J Appl Microbiol. 2009;106(1):34ā€“40.

    ArticleĀ  CASĀ  Google ScholarĀ 

  48. Sandegren L. Selection of antibiotic resistance at very low antibiotic concentrations. Ups J Med Sci. 2014;119:103ā€“7.

    ArticleĀ  Google ScholarĀ 

  49. Pereira R, Siler J, Ng J, Davis M, Grohn Y, Warnick L. Effect of on-farm use of antimicrobial drugs on resistance in fecal Escherichia coli of preweaned dairy calves. J Dairy Sci. 2014;97:7644ā€“54.

    ArticleĀ  CASĀ  Google ScholarĀ 

  50. Van Vleck PR, Lima S, Siler JD, Foditsch C, Warnick LD, Bicalho RC. Ingestion of Milk containing very low concentration of antimicrobials: longitudinal effect on fecal microbiota composition in Preweaned calves. PLoS One. 2016;11:e0147525.

    ArticleĀ  Google ScholarĀ 

  51. Thames CH, Pruden A, James RE, Ray PP, Knowlton KF. Excretion of antibiotic resistance genes by dairy calves fed milk replacers with varying doses of antibiotics. Front Microbiol. 2012;3:1ā€“12.

    ArticleĀ  Google ScholarĀ 

Download references

Acknowledgments

We thank Adam Roberts, Bastiaan Krom and Mark van Passel for critically reading the manuscript.

Funding

Financial support for this study was provided by the Stichting Kwaliteitsgarantie Vleeskalversector (SKV) and a precompetitive Dutch government EZ-Co-financing grant nr. 060.24960.

Availability of data and materials

The datasets supporting the conclusions of this article are available in the Sequence Read Archive (SRA) repository, project ID: PRJNA511412; https://www.ncbi.nlm.nih.gov/bioproject/PRJNA511412.

Author information

Authors and Affiliations

Authors

Contributions

BK, RM and VA assisted with interpretation of the data and wrote the manuscript; TvdB performed the statistical analyses of the qPCR data, and 16S data. MC assisted in the resistome analysis, MH performed the DNA isolation and metagenomic sequencing, MvanB and DM assisted in the antibiotic residue analysis, Eric Schoen conducted the statistical analysis of the resistome data, RB, AvdB, and LB were responsible for the intervention study and sample collection. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Bart J. F. Keijser or Tim J. van den Broek.

Ethics declarations

Ethics approval and consent to participate

All animal studies were performed in compliance with the Dutch law on animal experimentation (Wet op dierproeven, Wod / Directive 2010/63/EU of the European parliament and of the council of 22 September 2010 on the protection of animals used for scientific purposes). The study design was authorized by the Technical Committee of the SKV (Foundation for Quality Guarantee Veal Sector) on the 4th of February 2015, document TC15. Animals were placed under regular agricultural care according to the regulations of the Integrated Chain Management System (IKB) and regular clinical veterinary practices according to the general rules laid down in Council Directive 98/58/EC. Animal transport and handling during transport was compliant to EU directive 1/2005/EG and 1255/97/EG.

Consent for publication

All authors of the manuscript have read and agreed to its content and are accountable for all aspects of the accuracy and integrity of the manuscript in accordance with ICMJE criteria.

Competing interests

TNO is an independent non-profit research organization (employer of BK, TvdB, RM, VA, MC, MH, and ES) and has previously received grants from SKV. Ducares (employer of DM, MvB and AvdB) is a spin off company of TNO, and provides fee for service in biochemical analysis to both TNO and SKV. TNO holds a minority share in the company. Stichting Kwaliteitsgarantie Vleeskalversector (SKV) is an independent foundation, aiming to support for safe and responsibly produced veal.

Publisherā€™s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Figure S2. Spatial distribution of the pens where the calves were housed. The colors indicate the locations of the control group animals (red), the animals receiving a low dose of OTC (yellow), and the animals receiving a high dose of OTC (blue). (PNG 82 kb)

Additional file 2:

Figure S1. Experiment timeline. The upper part of the figure indicates the time points at which specific antibiotics were administered to the intervention and control groups. The lower part of the figure indicates the time points at which specific analyses were performed. (PNG 36 kb)

Additional file 3:

Table S1. Primers and probes used in qPCR. Table S2. A. Numbers of gene orthologues with significantly (P<0.01) different abundances between groups. Down and up indicate whether the genes are under- or over-represented respectively in one group compared to the other. B. Gene orthologues with significantly (P<0.01) different abundances between groups. Table lists for all significant gene orthologues, the KEGG reference and annotation, and the the Benjamini-Hochberg adjusted p-value (padj). Table S3. Resistance genes found in the metagenomic dataset that had a significantly higher abundance (p < 0.05) in one of the intervention groups compared to the control. (DOCX 26 kb)

Additional file 4:

Figure S3. Bar graph showing the taxonomic composition at the genus level in the control and intervention groups at different time points. The graph shows average abundance values for each group and time point. The ā€œrestā€ group includes genera detected with a relative abundance lower than 1%. (PNG 53 kb)

Additional file 5:

Figure S4. Heatmap of the ten operational taxonomic units identified through Canonical Analysis of Principal coordinates (CAP), (A) linked to calf age (time), (B) linked to group differences. (ZIP 20 kb)

Additional file 6:

Figure S5. Boxplots shoring the Shannon diversity of each group of calves at each time point. The boxes in the plot represent the interquartile ranges, the horizontal lines give the position of the medians, the vertical bars indicate the range. The dots indicate outliers. (PNG 329 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Keijser, B.J.F., Agamennone, V., van den Broek, T.J. et al. Dose-dependent impact of oxytetracycline on the veal calf microbiome and resistome. BMC Genomics 20, 65 (2019). https://doi.org/10.1186/s12864-018-5419-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-018-5419-x

Keywords