Animal sampling and follow-up
We collected fecal samples from 15 calves each from three veal farms in Brittany, France, during fattening. They were sampled by rectal swabbing at days 7 and 21 and then monthly for five months until the departure of the batch to the slaughterhouse. The calves were 14 days old when they arrived on the farms and were mainly fed milk replacer throughout the fattening period, which is reconstituted from milk powder with hot water. Calves stayed for 161 days on farms A and B and 147 days on farm C. Collective antibiotic treatments were recorded by the three farms throughout the fattening period (Fig. 1). Antibiotics were always used at therapeutic doses and administered orally in the feed. All calves received antibiotics more than once and calves from farms A and C received several consecutive antibiotic treatments during the first month (Fig. 1).
All calves from farms A and B included in the study were followed until the end of the fattening period. One calf from farm C died during fattening and was excluded from the study. For the 44 remaining calves, there were no missing samples during the follow-up and thus downstream analyses were performed on 308 samples.
After sampling, swabs were placed immediately in portable coolers with ice pack. Swabs were shipped to the antimicrobial resistance and virulence lab of the French Agency for Food, Environmental and Occupational Health & Safety (ANSES) lab in Lyon, France, within 24 hours, and stored at -80°C. After genomic DNA extraction from the swabs, we characterized the fecal microbiota by 16S rRNA gene amplicon sequencing (V4 region) using Illumina MiSeq technology. E. coli was quantified by qPCR targeting the 16S rRNA gene sequence specific to the Escherichia genus.
16S amplicon sequencing
After processing reads using the mothur pipeline, 34,153,188 quality amplicons were generated, with an average of 111,248 ± 63,002 per sample (Additional file 1: Fig. S1a). One sample had mostly poor-quality reads and a very small number of amplicons and was thus excluded from the downstream 16S sequencing analyses. The minimum number of operational taxonomic units (OTUs) detected in a sample was 119 and the maximum 1,302 (Additional file 1: Fig. S1b). The rarefaction threshold was set to 47,000 sequences (Additional file 1: Fig. S1c). Seven samples were below this threshold and excluded from the α- and β-diversity analyses.
Weighted and unweighted Unifrac distances
The similarity of microbiota composition among calves was tracked using β-diversity measures, which represents the dissimilarity between samples. The weighted Unifrac distances between calves were the highest on day 7, with an overall mean of 0.45 (± 0.10) (Fig. 2a). The weighted Unifrac distances started to decrease at the next sample, the overall mean reaching 0.33 (± 0.06) on day 21 (Fig. 2b). Weighted Unifrac distances remained low until the end of fattening, with a mean of 0.31 (± 0.08) on the last day (Fig. 2c). The time of sampling had a significant effect on the calf microbiome composition (p = 0.001, PERMANOVA on weighted Unifrac distances). It explained 15.5% of the between-sample variation in calves, indicating extensive sharing of the microbial community among calves at a given time of the fattening period. The microbiota of individuals belonging to the same farm were no more similar than those of individuals from different farms (p = 0.4, PERMANOVA on weighted Unifrac distances). Despite the use of different antibiotic molecules and different times of collective antibiotic treatments between farms, the farm affiliation explained only 4.0% of the between sample variation (Fig. 2).
We obtained similar results for the unweighted Unifrac distances (Additional file 2: Fig. S2, p = 0.001 and p = 0.5 for the time of sampling and farm in the PERMANOVA test, respectively). The farm affiliation explained 2.9% of the between-sample variation.
The weighted Unifrac distances between consecutive samplings showed a trend over time, with a downward trajectory for almost all calves (Additional file 3: Fig. S3). The mean intra-calf weighted Unifrac distances between the first and second samplings and the second and third samplings was 0.40 (± 0.10) and 0.33 (± 0.08), respectively. The mean intra-calf weighted Unifrac distances between the second to last sampling and last sampling was 0.25 (± 0.07). These results suggest that the magnitude of changes tended to decrease on a monthly scale.
Taxonomic composition of the microbiota
We investigated the taxonomic composition at different levels, from phyla to OTUs. We detected 19 phyla in the samples. Among them, only Firmicutes, Bacteroidetes, Actinobacteria, and Proteobacteria were present at a relative abundance above 1% in all samples. The phyla Firmicutes and Bacteroidetes were dominant on all farms throughout the period of the study. Their overall mean relative abundances were 47.7% (± 10.2) and 40.3% (± 11.3), respectively. The overall mean relative abundance of Actinobacteria was 5.2% (± 7.2) and that of Proteobacteria 3.7% (± 4.3). On farm A, Firmicutes was predominant at the beginning and then decreased slightly towards relative abundances similar to those of Bacteroidetes (Fig. 3).
We detected 349 genera, of which 50 (14%) were found to be one of the five most abundant taxa in at least one sample (Fig. 4, Additional file 4: Fig. S4). On day 7, the mean cumulative relative abundance of the five most abundant taxa was 73.9% (± 12.3) (Fig. 4a). On day 21, the mean cumulative relative abundance of the five most abundant taxa was 59.0% (± 7.4) and it was 60.2% (± 8.4) on the last day (Fig. 4b and Fig. 4c).
From the second month, 16 taxa at the genus level were detected in all calves at each sampling until the end of fattening. They were Alloprevotella, Bacteroides, Bifidobacterium, Blautia, Dorea, Faecalibacterium, Lactobacillus Olsenella, Parabacteroides, Prevotella, Pseudoflavonifractor,Ruminococcus2, and unclassified taxa from Clostridiales, Erysipelotrichaceae, Lachnospiraceae, and Ruminococcaceae. The overall mean cumulative relative abundance of these taxa was 67.5% (± 9.9), and was similar across samplings throughout the fattening period. On farm C, the mean was equal to 74.1% (± 5.5), 72.8% (± 5.5), and 70.0% (± 6.5) on days 35, 91 and 147, respectively. On farms A and B, the mean was 64.6% (± 10.0), 64.6% (± 8.8), and 65.5% (± 7.7) on days 49, 106, and 161, respectively.
For each calf, the proportion of OTUs not previously detected was higher than the proportion of OTUs detected in the previous sample (Fig. 5). The proportion of OTUs detected in the previous sample varied between 7.2% and 48.6%, whereas the proportion of OTUS that had never been detected before varied between 51.4% and 92.8%. These results suggest that the temporal dynamics of the calf microbiota is driven by the replacement of autochthonous OTUs by the newly detected ones. The proportion of newly detected OTUs tended to decrease over time, concomitantly with an increase in the proportion of OTUs detected in two consecutive samples (Fig. 5). Only 50 OTUs simultaneously persisted in more than 97% of calves between consecutive samples. Two OTUs from the genera Faecalibacterium and an unclassified OTU from the Ruminococcaceae family persisted from the first to last month in more than 97% of calves (Additional file 5: Table S1). No OTU was newly detected or lost simultaneously by more than 50% of the calves (Additional file 5: Table S1).
No over-abundance or depletion of taxa was detected at the phylum (Fig. 3), nor genus level (Fig. 4a and Fig. 4b) in samples of calves under antibiotic treatment or those that had received antibiotics in the 15 days before sampling relative to calves not exposed during the same period. For example, at day 21, although calves on farm A had received antibiotics for 20 days while calves on farm B had not, 154 genera were detected on both farms. They represented 86.5% and 85.6% of the genera detected that day on farms A and B, respectively. At day 21, 157 genera were detected on both farms B and C. They represented 87.2% and 86.7% of the genera detected that day on farms B and C, respectively, despite the fact that calves on farm C had received antibiotics for 14 days. Sixteen taxa at the genus level were only detected on farm B and were found at a relative abundance of < 1% (Chlamydophila, Snodgrassella, unclassified taxa from Rhodobacteraceae, Tissierella, Clostridium_XI, Comamonas, Basfia, Janibacter, Anaerosporobacter, Pseudoscardovia, unclassified taxa from Deltaproteobacteria, Brevinema, Brucella, unclassified taxa from Gammaproteobacteria, unclassified taxa from Desulfovibrionaceae, Oligosphaera).
Shannon index and the number of observed OTUs
We evaluated the temporal dynamics of microbiota diversity in each sample by examining the α-diversity metrics, the Shannon index, and the number of observed OTUs. The Shannon index showed an increasing trend over time, with an overall mean of 3.05 (± 0.82) on day 7 and 4.23 (± 0.59) at the end of fattening. The temporal dynamics of microbiota diversity were best described by a two-slope model (Likelihood Ratio Test between the two candidate models, p < 10-15, Fig. 6a, Additional file 6: Table S2). The coefficient of the first slope was higher than that of the second, suggesting a large increase in diversity during the first month, followed by a lowering of the rate of increase, and even stagnation on farm A. The estimates for the intercept and second slope for farms B and C were significantly higher than those for farm A. There was a similar temporal trend and a significant farm effect on the number of observed OTUs (Additional file 7: Fig. S5a, Additional file 6: Table S3).
A variable representing antibiotic treatment in the 15 days previous to sampling (still ongoing or not at the time of sampling) was added to the two-slope model and tested for significance for both indices. There was a significant effect of antibiotics on the Shannon index intercept, with an estimated decrease of -0.17 CI95%[-0.27; -0.06] (p = 0.003, Additional file 6: Table S2, Fig. 6b), indicating an antibiotic-induced reduction of bacterial diversity during exposure. The effect of antibiotic exposure on the number of observed OTUs was also significant and showed a similar effect of decreasing diversity, with an estimated decrease of 82.7 OTUs CI95%[-115.8; -49.6] (p = 1*10-6, Additional file 6: Table S3, Additional file 7: Fig. S5b).
Absolute number of E. coli per gram of feces
We quantified the absolute number of Escherichia/g by qPCR, which can be considered as a fair proxy of the absolute number of E. coli in calves (see Methods). E. coli is the main facultative anaerobic bacteria in the large intestine and a marker of dysbiosis [12]. We first compared these absolute numbers to the relative abundance of the Escherichia genus estimated by 16S rRNA gene sequencing. The number of E. coli estimated by qPCR was strongly and positively associated with the relative abundance of the Escherichia genus (Spearman’s correlation, rs = 0.80, p < 10-15, Fig. 7a), confirming the strong relationship between the two variables. The temporal dynamics of the number of E. coli/g was similar for the three farms, with an overall mean of 7.81 log10 (E. coli/g) (± 0.67) on day 7. During the second month, a transient but important increase (approximately 2 log10) occurred on the three farms (Fig. 7b). The dynamics of the number of E. coli/g are best described by a quartic function of time (LRT between cubic and quartic function of time models, p = 0.01). The final model gave the same estimates for the parameters of farms A and C, suggesting the presence of an additional factor in shaping the number of E. coli/g on these farms (Fig. 7b and Additional file 6: Table S4).
As for α-diversity indices, a variable representing antibiotic treatment in the 15 days prior to sampling (still ongoing or not at the time of sampling) was added to the quartic model and tested for significance. There was a significant reduction in the E. coli population in calves treated in the previous 15 days relative to those that were not, with an estimated decrease of -0.37 log10 (E. coli/g) CI95%[-0.66; -0.08] (p = 0.01, Additional file 6: Table S4, Fig. 7c).
Association between the estimated dose of milk powder given in two farms and bacterial abundance
In farms B and C, for which the information was available, we explored associations between the relative abundance of genera and the daily dose of milk powder recommended by the integrator. Calves were almost exclusively fed milk replacer throughout the fattening period, which is reconstituted from milk powder with hot water. Their diet was also supplemented with a small amount of solid feed from the first weeks. We conducted 349 Spearman rank correlation tests and found 17 taxa at the genus level for which the Spearman correlation coefficient was positive and significantly different from zero (Table 1, Additional file 8: Table S5). The genus with the strongest correlation was Megasphaera, with a moderate to strong association (rs = 0.60, Bonferroni adjusted p < 1*10-15, Table 1, Additional file 9: Fig. S6). The genera with the next highest positive correlation coefficients were Enterococcus, Dialister, and Mitsuokella, indicating moderate association with the dose of milk powder (rs = 0.44, 0.42, and 0.41, and Bonferroni adjusted p = 2*10-8, 3*10-7, and 6*10-7, respectively, Table 1, Additional file 9: Fig. S6). The genus Escherichia was also positively associated with milk powder, but to a lesser extent (rs = 0.28, Bonferroni adjusted p = 0.02, Table 1).
Table 1. Correlation between the dose of milk powder and fecal microbiota on farms B and C.
Taxa identified at the genus level*
|
Spearman correlation coefficient**
|
p-value
|
Bonferroni adjusted p-value
|
Megasphaera
|
0.60
|
2*10-20
|
6*10-18
|
Enterococcus
|
0.44
|
5*10-11
|
2*10-8
|
Dialister
|
0.42
|
1*10-9
|
3*10-7
|
Mitsuokella
|
0.41
|
2*10-9
|
6*10-7
|
Collinsella
|
0.37
|
9*10-8
|
3*10-5
|
Unclassified Veillonellaceae
|
0.32
|
4*10-6
|
0.001
|
Sutterella
|
0.31
|
1*10-5
|
0.003
|
Selenomonas
|
0.31
|
1*10-5
|
0.003
|
Anaerofilum
|
0.30
|
2*10-5
|
0.006
|
Faecalibacterium
|
0.30
|
2*10-5
|
0.007
|
Butyricimonas
|
0.29
|
4*10-5
|
0.01
|
Megamonas
|
0.28
|
6*10-5
|
0.02
|
Escherichia/Shigella
|
0.28
|
8*10-5
|
0.02
|
Fastidiosipila
|
0.27
|
0.00013
|
0.038
|
Helcococcus
|
0.27
|
0.00014
|
0.040
|
Bifidobacterium
|
0.27
|
0.00014
|
0.041
|
Psychrobacter
|
0.27
|
0.00015
|
0.045
|
* Data represent the taxa identified at the genus level with a positive correlation with the dose of milk powder estimated on farms B and C. Only taxa with a significant correlation (Bonferroni corrected p-value < 0.05) are shown.
** Spearman correlation coefficients ∈ [0.4; 0.6] indicate moderate to strong associations, whereas coefficients < 0.4 indicate weak associations and were not considered for further analyses.
As E. coli is a β-galactosidase-positive species, we compared the estimated dynamic profile of the absolute number of E. coli/g and the dose of milk powder recommended by the integrators on farms B and C. The estimated profiles were very close for both farms, particularly during the second month, with the peak of the dose of milk powder superimposed over that for the number of E. coli/g (Fig. 7d). We searched for an association between the estimated daily doses of milk powder and the farm predicted numbers of E. coli/g on farms B and C using Spearman’s correlation test. We found a significant strong positive association between the farm predicted numbers of E. coli/g and the estimated dose of milk powder (rs = 0.77, p< 1*10-15).