# Temporal dynamics of the fecal microbiota in veal calves in a 6-month field trial

Méril Massot, Marisa Haenni, Thu Thuy Nguyen, Jean-Yves Madec, France

Mentré, Erick Denamur

► To cite this version:

Méril Massot, Marisa Haenni, Thu Thuy Nguyen, Jean-Yves Madec, France Mentré, et al. Temporal dynamics of the fecal microbiota in veal calves in a 6-month field trial. Animal Microbiome, 2020, 2 (1), pp.32. 10.1186/s42523-020-00052-6. anses-03992681

HAL Id: anses-03992681

https://anses.hal.science/anses-03992681

Submitted on 16 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Temporal dynamics of the fecal microbiota in veal calves in a 6-month field trial

<figure><img src="image_2.png" /></figure>

Check for
updates

Méril Massot¹*, Marisa Haenni², Thu Thuy Nguyen¹, Jean-Yves Madec², France Mentré¹,³ and Erick Denamur¹,⁴

## Abstract

**Background:** Little is known about maturation of calves' gut microbiome in veal farms, in which animals are confined under intensive-farming conditions and the administration of collective antibiotic treatment in feed is common. We conducted a field study on 45 calves starting seven days after their arrival in three veal farms. We collected monthly fecal samples over six months and performed 16S rRNA gene sequencing and quantitative PCR of *Escherichia coli* to follow the dynamics of their microbiota, including that of their commensal *E. coli* populations. We used mixed-effect models to characterize the dynamics of α-diversity indices and numbers of *E. coli*, and searched for an effect of collective antibiotic treatments on the estimated parameters. On two farms, we also searched for associations between recommended daily doses of milk powder and bacterial abundance.

**Results:** There was high heterogeneity between calves' microbiota upon their arrival at the farms, followed by an increase in similarity, starting at the first month. From the second month, 16 genera were detected at each sampling in all calves, representing 67.5% (± 9.9) of their microbiota. Shannon diversity index showed a two-phase increase, an inflection occurring at the end of the first month. Calves receiving antibiotics had a lower intercept estimate for Shannon index (− 0.17 CI<sub>95%</sub>[-0.27; − 0.06], *p* = 0.003) and a smaller number of *E. coli*/ gram of feces during the treatment and in the 15 days following it (− 0.37 log<sub>10</sub> (*E. coli*/g) CI<sub>95%</sub>[-0.66; − 0.08], *p* = 0.01) than unexposed calves. There were moderate to strong positive associations between the dose of milk powder and the relative abundances of the genera *Megasphaera*, *Enterococcus*, *Dialister* and *Mitsuokella*, and the number of *E. coli* (*r*<sub>s</sub> ≥ 0.40; Bonferroni corrected *p* < 0.05).

**Conclusions:** This observational study shows early convergence of the developing microbiota between veal calves and associations between the dose of milk powder and members of their microbiota. It suggests that administration of collective antibiotic treatment results in a reduction of microbial diversity and size of the *E. coli* population and highlights the need for additional work to fully understand the impact of antibiotic treatment in the veal industry.

**Keywords:** Veal calves, Fecal microbiota development, *Escherichia coli*, 16S rRNA gene sequencing, Milk powder, Antibiotics

## Background

The veal-calf industry is an intensive farming system that produces meat from milk-fed calves. Europe is one of the largest producer of veal in the world, producing approximately six million heads yearly, mostly in France,

the Netherlands, Italy, and Belgium. Male dairy calves are commonly used to produce veal. They are purchased in dairy farms by “integrators”, companies that are involved in all stages of the production process [1]. They are collected at 2 weeks of age, batched, and placed without delay in dedicated closed buildings for approximately 6 months. Calves are mainly fed milk replacers and a small amount of solid feed is introduced during the first weeks for the welfare of the animals. They are

<footnote>* Correspondence: meril.massot@inserm.fr
¹Université de Paris, IAME, INSERM, Site Xavier Bichat, 16 rue Henri Huchard, F-75018 Paris, France
Full list of author information is available at the end of the article</footnote>
reared indoors, often in high-animal-density buildings, which increases the risk of pathogens spreading among the batch. Collective antibiotic treatment is frequent, particularly at the start of the fattening period, when calves coming from different dairy farms are grouped together. In France, veal calves typically receive more than eight antibiotic treatments during the fattening process, mainly by the oral route [2, 3]. Collective treatment represents almost 96% of the administrated treatments and 88% of the prescriptions are made in advance of an anticipated outbreak of disease, when some animals in the batch develop symptoms of infection [2–4]. A pervasive effect of antibiotic treatment can be the collapse of gut bacterial populations, which results in a loss of fecal microbial diversity, as shown in cattle [5, 6].

Because of the extensive use of antibiotics, combined with specific dietary practices and intensive farming practices, maintaining a proper gut microbiota balance during calf development is a key challenge in the veal-calf industry. To our knowledge, only a few attempts have been made to investigate the effect of antibiotics on the gut microbiota composition of pre-weaned dairy calves in commercial farms [7, 8], and no study has focused on calves reared in veal farms.

Based on a field trial, we report the temporal dynamics of bacterial communities in veal calves highly exposed to antibiotics from an early age, with an additional focus on the commensal *E. coli* population. Although *E. coli* is one of the most common pathogens that causes diarrhea in calves, it is also found as a commensal in the gut of healthy calves [9, 10]. Few studies have focused on commensal *E. coli* populations in pre-weaned calves [11]. We hypothesized that the developmental trajectory of the microbiota is influenced by the use of antibiotics during growth, resulting in the development of a community with lower diversity and persistent shifts in taxonomic composition. We addressed this question by performing 16S rRNA gene sequencing and quantitative PCR (qPCR) of *Escherichia* on 312 rectal swabs collected from 45 veal calves distributed in French veal farms over 6 months.

# Results

## Animal sampling and follow-up

We collected fecal samples from 15 calves each from three veal farms partnered with different integrators in Brittany, France, during fattening. They were sampled by rectal swabbing at days 7 and 21 and then monthly for 5 months until the departure of the batch to the slaughter-house. 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 packs. 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 h, 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 amplicon sequencing analyses. The minimum number of operational taxonomic units (OTUs) detected in a sample was 119 and the maximum 1302 (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, permutational multivariate analysis of variance (PERM ANOVA) on weighted Unifrac distances). It explained
<figure><img src="image_3.png" /></figure>

**Fig. 1** Scheme of sampling dates for each farm. Sampling points for farm A, farm B, and farm C are represented in the upper panel, middle panel, and lower panel, respectively. “N” indicates the number of calves studied on each farm. The days of sampling are indicated by grey dots. Of note, one calf on farm C died during the fattening period and was excluded from the study. Antibiotic treatments are indicated by bold dark lines or back triangles, and the names of the antibiotics are given in the legend

<figure><img src="image_4.png" /></figure>

**Fig. 2** Heatmaps of the weighted Unifrac distances at the first month and last month of fattening. Heatmaps of the β-diversity weighted Unifrac distances matrix are shown for the (a) first sampling (day 7), (b) second sampling (day 21), and (c) last sampling (day 161 for farms A and B and day 147 for farm C). Each square represents a pairwise distance between two calves. The pale yellow squares indicate low Unifrac distances, whereas dark red squares indicate high Unifrac distances. The calves are ordered according to the farms in both the lines and columns. The calf's distance from itself is represented by the white square on the main diagonal. The means ± standard deviations for each sampling and farm are shown in the lower triangles
15.5% of the between-sample variation in calves, indicat-
ing extensive sharing of the microbial community among
calves at a given time of the fattening period. The micro-
biota 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 dis-
tances). Despite the use of different antibiotic molecules
and different times of collective antibiotic treatments be-
tween farms, the farm affiliation explained only 4.0% of
the between sample variation (Fig. 2).

We obtained similar results for the unweighted Uni-
frac distances (Additional file 2: Fig. S2, p = 0.001 and
p = 0.5 for the time of sampling and farm in the PERM
ANOVA test, respectively). The farm affiliation ex-
plained 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), respect-
ively. The mean intra-calf weighted Unifrac distances be-
tween 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 c).

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*, 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% (±

<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /></figure>

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. An OTU from the genus *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).
<figure><img src="image_7.png" /></figure>

No over-abundance or depletion of taxa was detected at the phylum (Fig. 3), nor genus level (Fig. 4a and b) 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).

<figure><img src="image_8.png" /></figure>
## 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 (LRT) between the two candidate models, *p* < 10⁻¹⁵, 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 CI<sub>95%</sub>[-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 CI<sub>95%</sub>[-115.8; -49.6] (*p* = 1*10<sup>-6</sup>, Additional file 6: Table S3, Additional file 7: Fig. S5b).

<figure><img src="image_9.png" /><figcaption>**Fig. 7** Absolute quantification of the *Escherichia coli* population by qPCR. **a** Relative abundance of the *Escherichia* genus as a function of the number of *E. coli/g* estimated by qPCR. Each point represents a sample. Points on the x-axis represent samples for which no 16S rRNA gene sequence of the *Escherichia* genus was detected by sequencing. **b** Dynamics of the mean observed and predicted number of *E. coli/g* for each farm in the final model without the antibiotic-treatment effect. The mean values ± standard deviations of the observed data for each farm at each sampling time are represented by the dashed bars. The model-predicted profiles and their 95% confidence bands are represented by the solid lines and bands, respectively. The predicted profiles of farms A and C are overlapping, which is why the predicted profile of farm A does not appear. **c** Dynamics of the mean observed and predicted number of *E. coli/g* for each farm in the final model with the antibiotic-treatment effect. Antibiotic treatments during sampling or within 15 days before sampling are color-coded by farm and indicated above the x-axis. **d** Temporal dynamics of the recommended dose of milk powder per kilo of live weight for farms B and C</figcaption></figure>
## 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' feces (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, r<sub>s</sub> = 0.80, *p* < 10<sup>-15</sup>, 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 log<sub>10</sub> (*E. coli/g*) (± 0.67) on day 7. During the second month, a transient but important increase (approximately 2 log<sub>10</sub>) occurred on the three farms (Fig. 7b). The dynamics of the number of *E. coli/g* were 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 log<sub>10</sub> (*E. coli/g*) CI<sub>95%</sub>[-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 (r<sub>s</sub> = 0.60, Bonferroni adjusted *p* < 1*10<sup>-15</sup>, 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

<table><caption>**Table 1** Correlation between the dose of milk powder and fecal microbiota on farms B and C</caption><thead><tr><th>Taxa identified at the genus level*</th><th>Spearman correlation coefficient<sup>a</sup></th><th>p-value</th><th>Bonferroni adjusted p-value</th></tr></thead><tbody><tr><td>Megasphaera</td><td>0.60</td><td>2*10<sup>-20</sup></td><td>6*10<sup>-18</sup></td></tr><tr><td>Enterococcus</td><td>0.44</td><td>5*10<sup>-11</sup></td><td>2*10<sup>-8</sup></td></tr><tr><td>Dialister</td><td>0.42</td><td>1*10<sup>-9</sup></td><td>3*10<sup>-7</sup></td></tr><tr><td>Mitsuokella</td><td>0.41</td><td>2*10<sup>-9</sup></td><td>6*10<sup>-7</sup></td></tr><tr><td>Collinsella</td><td>0.37</td><td>9*10<sup>-8</sup></td><td>3*10<sup>-5</sup></td></tr><tr><td>Unclassified Veillonellaceae</td><td>0.32</td><td>4*10<sup>-6</sup></td><td>0.001</td></tr><tr><td>Sutterella</td><td>0.31</td><td>1*10<sup>-5</sup></td><td>0.003</td></tr><tr><td>Selenomonas</td><td>0.31</td><td>1*10<sup>-5</sup></td><td>0.003</td></tr><tr><td>Anaerofilum</td><td>0.30</td><td>2*10<sup>-5</sup></td><td>0.006</td></tr><tr><td>Faecalibacterium</td><td>0.30</td><td>2*10<sup>-5</sup></td><td>0.007</td></tr><tr><td>Butyricimonas</td><td>0.29</td><td>4*10<sup>-5</sup></td><td>0.01</td></tr><tr><td>Megamonas</td><td>0.28</td><td>6*10<sup>-5</sup></td><td>0.02</td></tr><tr><td>Escherichia/Shigella</td><td>0.28</td><td>8*10<sup>-5</sup></td><td>0.02</td></tr><tr><td>Fastidiosipila</td><td>0.27</td><td>0.00013</td><td>0.038</td></tr><tr><td>Helcococcus</td><td>0.27</td><td>0.00014</td><td>0.040</td></tr><tr><td>Bifidobacterium</td><td>0.27</td><td>0.00014</td><td>0.041</td></tr><tr><td>Psychrobacter</td><td>0.27</td><td>0.00015</td><td>0.045</td></tr></tbody></table>

<footnote>* 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</footnote>

<footnote><sup>a</sup> 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</footnote>
of milk powder ($r_s = 0.44, 0.42, \text{ and } 0.41$, and Bonferroni adjusted $p = 2*10^{-8}, 3*10^{-7}, \text{ 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 ($r_s = 0.28$, Bonferroni adjusted $p = 0.02$, Table 1).

As *E. coli* is a β-galactosidase-positive species, which means it is able to cleave lactose into monosaccharides, we compared the estimated dynamic profiles 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 ($r_s = 0.77, p < 1*10^{-15}$).

## Discussion

We characterized the dynamics of the fecal microbiota of calves from two weeks to six months of age on three commercial veal farms representative of the three main French integrators and of management practices in the veal industry in France. The calves were mainly fed milk replacers throughout the follow-up and received several collective antibiotic treatments at therapeutic doses, most administered in the first weeks of fattening (Fig. 1). We performed 16S rRNA gene sequencing to study the composition of the microbiota and qPCR of the *Escherichia* genus as a proxy of *E. coli* to quantify its commensal populations. We estimated the daily dose of milk powder recommended by the integrator for two farms to search for an association with the relative abundance of the detected genera. The most striking results of this study are (i) the convergence of the fecal microbiota composition among calves, which began during the first month of life, along with an increase in α-diversity, (ii) a decrease in microbiota diversity and the size of the *E. coli* population during or within the 15 days following an antibiotic treatment relative to non-exposed calves of the same age (reduction of the Shannon index by 0.17 and the number of *E. coli/g* of feces by 0.37 $log_{10} (E. coli/g)$), and (iii) a significant association between the estimated daily dose of milk powder and the relative abundance of four genera (*Megasphaera*, *Enterococcus*, *Dialister*, *Mitsuokella*) and the predicted farm profiles of the number of *E. coli/g* from our model.

The development of the microbiota of these calves was characterized by the dominance of a small number of genera mainly from the *Firmicutes* and *Bacteroidetes*

phyla (Fig. 3). These phyla remained dominant as the calves aged (Fig. 4), with a simultaneous increase in microbiota diversity (Fig. 6 and Additional file 7: Fig. S5a). These developmental features have already been described for calves with the same characteristics (age, sex, and breed) fed milk replacers [10], as well as for females of the same breed in Canada [10, 13], the USA [14, 15], and Japan [16] and for dairy calves of a different breed in Austria [17]. These shared findings suggest that the fecal microbiota of calves undergoes a predictable age-dependent trend that is common to distinct calf populations.

The high heterogeneity of the microbiota composition at day 7 on the farms (corresponding to 3 weeks of age) is likely attributable to the distinct origin of the calves, as they came from different dairy farms. Prior studies have noted the important influence of exposure to the bacterial communities of both the dam and the environment on the composition of the microbiota throughout the gut of the newly born calf [18, 19]. The transport from dairy farms to veal farms could also be responsible for the high heterogeneity between calves at the beginning of fattening. Transport has been reported to have disrupting effect on the gut microbiota in young beef cattle 5 days after their transport to the feedlot [6]. The composition of milk replacers may have been a source of variability between farms during fattening. A difference in the relative abundance of *Bifidobacterium* spp. and *Faecalibacterium prausnitzii* was found between one-month-old dairy calves fed milk replacers containing different levels of proteins and fat, although this difference was transitory [20].

From as soon as the end of the first month after arriving on the farms, bacterial succession (Fig. 5) gradually increased the similarity of the microbiota composition among the calves on all farms. This increase was due to both an increase in sharing of the same bacterial members and the homogeneity of their relative abundance, as shown by unweighted and weighted Unifrac distances (Additional file 2: Fig. S2 and Fig. 2, respectively). This convergent pattern occurred in the absence of environmental or dietary changes, such as weaning, the calves being reared in dedicated closed buildings and drinking milk replacers throughout fattening. This suggests that the influence of environmental and dietary factors on such convergence was probably limited and highlights the likely role of host physiology. Convergence related to age has also been observed in the ruminal microbiota of calves between one day and two years of age [21] and in both the ruminal and fecal microbiota of dairy heifers receiving different diets before weaning [22]. These studies suggest that such changes in composition are not restricted to the lower part of the gut and are not strongly driven by diet. An intriguing possibility is that such
convergent stabilization of the microbiota composition over time (Fig. 2, Additional file 2: Fig. S2) may be linked to age-dependent shifts of the gut mucosal immune system, as the expression of Toll-like receptors in both the rumen and colon have been shown to change as calves age [23]. Constraints imposed by the gut environment and autochthonous microbiota on allochthonous bacterial settlement may become less permissive, resulting in more specific requirements as the calves aged.

The microbiota of calves on farms where collective antibiotic treatments were given in the previous 15 days or during sampling underwent a reduction in diversity and the number of *E. coli* relative to calves of the same age that had not been exposed during the same period (Fig. 6b, Additional file 7: Fig. S5, Fig. 7c). We pooled the effects of both long-term and short-term antibiotic treatment and molecules with different spectra to focus on the common disruptive effects of antibiotics on microbial ecosystems. As all calves of the same farm received the same antibiotic treatment at the same time, the design of our study did not allow the analysis of the specific effects of each molecule. As the treatments were collective, we did not have negative controls for calves receiving antibiotics from the same farm. Hence, the results could have been biased by farm-specific variables, which could have modulated the relative abundance of certain genera. To our knowledge, no study has yet explored the impact of antibiotic treatments on the maturation of gut microbiota in veal calves. Moreover, literature on the effect of each antibiotic class on the hindgut microbiota of preweaned dairy calves is scarce. As far as we know, research on this topic has focused on a few classes: macrolides [7, 24], tetracycline [7], amphenicol [7], fluoroquinolones [24], and polypeptides [25]. These studies, with the exception of that by Xie et al. [25], all focused on the effects of a single injection of the antibiotic, which differs from the treatment plans and route of administration (oral) in our study. Xie et al. reported dysbiosis in neonatal dairy calves receiving bacitracin for 10 days in the milk replacer, characterized by an increase in the abundance of *Escherichia* and *Enterococcus*, along with a decrease in the abundance of *Dorea*, *Collinsella*, *Eubacterium*, *Faecalibacterium*, *Papillibacter*, *Peptostreptococcus*, *Prevotella*, and *Roseburia*, relative to non-exposed calves of the same age [25]. The modalities of treatment in this study were close to those for the calves of farms A and C just after their arrival on farms. On farm A, they received a combination of antibiotics in their feed for several days, including colistin, another polypeptide. Their finding concerning the genus *Escherichia* is contrary to the negative effect of antibiotics that we found. This may be explained by the difference of the spectrum between bacitracin and the antibiotics used on the veal farms, as bacitracin has a

narrow spectrum, targeting Gram-positive bacteria, whereas *E. coli* is sensitive to antibiotics used during extended treatment (colistin, sulfonamides, tetracycline, trimethoprim). It may also be partially explained by the age of calves, as the microbiota change markedly between birth and the first month of life [16, 26].

Antibiotic-induced loss of diversity has already been reported in young beef cattle [6], as well as pre-weaned calves [7], and was often found to be associated with the depletion of beneficial bacteria and/or the increase of opportunistic pathogens [7, 25]. Antibiotic-induced dysbiosis is also observed when pre-weaned calves are fed low doses of antibiotic molecules [27].

The effect of the antibiotic treatments was small relative to the longitudinal changes (Additional file 6: Tables S2, S3 and S4), consistent with findings in pre-weaned dairy calves receiving enrofloxacin and tulathromycin metaphylactic treatment [24] and beef cattle receiving oxytetracycline or tulathromycin injection [6]. These findings can be explained by the existence of a natural resistome, independently of any antibiotic treatment, carried by certain abundant families in the fecal microbiota of pre-weaned calves, as recently shown [28]. Generation of these antibiotic-resistance gene carrying families were found to be dominant in the feces of the veal calves, such as *Anaerostipes*, *Blautia*, and *Roseburia* (*Lachnospiraceae* family), *Enterococcus* (*Enterococcaceae* family), *Faecalibacterium* and *Pseudoflavonifractor* (*Ruminococcaceae* family), *Bacteroides* (Bacteroidaceae family), and *Streptococcus* (*Streptococcaceae*). Members of the *Enterobacteriaceae* family, such as *E. coli*, were also found to be a major reservoir of antibiotic-resistance genes within the microbiota resistome.

Two recent studies reported that the resistome of fecal microbiota in pre-weaned dairy calves is composed of resistance-conferring genes against tetracycline, sulfonamides, trimethoprim, β-lactams, and macrolides [8, 28]. These results suggest the existence of a natural resilience in fecal bacterial communities to collective antibiotic treatments on veal farms; the antibiotics used to treat the calves in our study belonged to these classes. None of these studies reported the presence of the *mcr-1* gene, which confers resistance to colistin, another molecule used to treat calves, in microbial communities, although another study detected the gene in the commensal *E. coli* of veal calves [3]. Moreover, it has been shown that this natural resistome is shaped by the bacterial phylogeny of the fecal microbiome and decreases as the calves age. One of the main drivers of such a decrease was the decrease in abundance of the *Enterobacteriaceae* family, in which 90% of the members were classified as *E. coli* [28]. It is well-known that commensal *E. coli* populations of veal calves harbor high levels of antibiotic-resistance genes [3, 29], and that they are diverse.
Antibiotic treatment may have promoted an increase in the number of specific pre-existing *E. coli* strains at the beginning of fattening, as the extended treatments during the first month did not result in marked depletion of the *E. coli* population.

Among others, we found a possible link between *E. coli* population dynamics and the use of milk replacer, which is reconstituted from dry milk powder and rich in lactose (Fig. 7d) [30]. Lactose allows the growth of the vast majority of *E. coli* strains [31]. It has been shown that the lag time and generation time of *E. coli* strains, which depend on metabolic efficiency and are crucial for gut colonization and persistence [32, 33], are influenced by the type and abundance of available nutrients in the habitat [34]. Other components of milk replacers may have also influenced the *E. coli* population dynamics, such as vitamin D, for which the absence of the vitamin D receptor in the intestinal epithelium of mice has been associated with increased *E. coli* loads [35]. Our findings are consistent with those of another study in which the fecal microbiota of Simmental calves was followed during their first 3 months of life, although only six calves were included in the study and the sampling was sparse [17]. The relative abundance of the genus *Escherichia* was found to be maximal during the milk-feeding period and decreased before weaning.

Lactose, which is one of the main components of milk powder (approximately 45% of the dry weight), is hydrolyzed to the monosaccharides glucose and galactose by bacteria that synthesize the β-galactosidase enzyme. We therefore looked for the presence of this enzyme in available annotated genomes from the NCBI genome database of the *Megasphaera*, *Enterococcus*, *Dialister*, and *Mitsuokella* genera. The sequence of *E. coli* LacZ β-galactosidase was also compared to the protein sequences found in members of these genera using the blastp program [36, 37]. A β-galactosidase gene sequence was found in the genome of *Mitsuokella multacida*, which was isolated from human feces [38], and several species of the *Enterococcus* genus, which have been isolated from cattle (*E. faecalis*, *E. faecium*, *E. hirae*, *E. thailandicus*, *E. malodoratus*, *E. devriesei*, *E. casseliflavus*, *E. italicus*) [39–42] (data not shown). The concomitant fluctuations of this lactose-rich source and the relative abundance of the genera *Enterococcus* and *Mitsuokella* strongly suggest a direct role of host diet on members of the fecal microbiota. We did not find a β-galactosidase annotated gene in the 24 *Megasphaera* annotated genomes nor in the 18 *Dialister* annotated genomes available on NCBI [37]. The bacteria of these two genera did not carry a protein similar to the *E. coli* LacZ β-galactosidase (data not shown). The discrepancy between the associations found for the genera *Megasphaera* and *Dialister* and the absence of the β-

galactosidase enzyme sequence in the genome of known members of these genera could be explained by the utilization of another nutrient present in the milk powder by the members of these taxa. Another explanation for this discrepancy could be the limited redundancy of carbon source use among members of these genera coupled with a small number of genomes from these genera in the NCBI database.

Our study had several limitations. First, we followed the fecal microbiota of calves reared in commercial veal farms. Thus, the calves were not randomly assigned to the different farms (nor antibiotic treatments) and neither the environmental nor dietary variables could be controlled as they would be in a randomized trial. Nevertheless, one of the aims of this study was to characterize the fecal microbiota of calves reared following common veal-farm practices. These three farms, in which field studies had already been conducted [3], are representative of management practices in the veal-calf industry in France. Second, the calves were only sampled after spending 7 days on the farms. Thus, no information in terms of their microbiota composition before antibiotic treatment was available. They were also sampled on a monthly basis, whereas high-frequency sampling has been recommended in early-life microbiome studies in infants [43]. Furthermore, sampling was performed independently of antibiotic treatment. Hence, certain short- and mid-term age- and antibiotic-associated changes may have been missed. Third, although we tracked the dynamics of microbiota using 16S rRNA gene sequencing and *E. coli* qPCR, we only focused on specific features of this complex ecosystem and may have missed specific patterns at other levels. For example, we had no information concerning the dynamics of fecal bacterial loads, which have been shown to vary in newborn calves [9, 44, 45] and be associated with microbiota composition [46].

Nevertheless, veal calves, as studied here, have attributes relevant to the exploration of the microbiota maturation process. First, batches are composed of male calves of the same age and breed (usually Holsteins). Thus, they are highly genetically and physiologically homogeneous. Second, they share the same living environment and diet, which are not subject to major changes, as they are reared in dedicated closed buildings in which the conditions are stable and controlled to optimize their growth. Moreover, they do not experience any drastic changes in their diet, as it remains predominantly composed of milk replacers during the 6 months of fattening. Third, the systematic administration of antibiotics at therapeutic doses to all members of the batch is common practice to prevent the spread of infectious diseases [2]. As healthy young subjects that share similar controlled conditions over a long period of time and experience common antibiotic exposure, veal calves represent a unique opportunity to disentangle the
factors that drive microbiota assembly under real-life conditions.

## Conclusion

This observational study conducted on calves reared under intensive-farming conditions shows (i) early convergence of the developing fecal microbiota among farms and (ii) a significant association between the estimated daily doses of milk powder and the relative abundance of certain genera and the predicted farm profiles of the number of *E. coli/g*. This study also suggests that the administration of collective antibiotic treatment results in a limited reduction of diversity and size of the *E. coli* population and highlights the need for additional studies to fully understand the impact of antibiotic treatment in the context of the veal industry. To our knowledge, this is the first field study to follow the microbiota composition and size of the commensal *E. coli* population of veal calves throughout the fattening period.

## Methods

### Animal handling and sampling

We collected fecal samples from veal calves during a cohort study dedicated to monitor the excretion of extended spectrum β-lactamase- (ESBL)-producing *E. coli*, which has been shown to be frequent in veal calves [3, 47]. The fecal excretion of ESBL-producing *E. coli* was followed in 45 veal calves distributed in three French veal farms (named A, B, and C) located in the region of Brittany [3], within a 100-km radius around Rennes. The three farmers raised calves in partnership with different integrators, which were the main veal calves producers in France. We streaked swabs on selective ChromID ESBL agar (bioMérieux, Marcy l'Etoile, France) and classified calves as "ESBL-producing *E. coli* high-level excretor", "low level excretor" or "ESBL-producing *E. coli*-free" based on the number of colonies that grew after 24 h at 37 °C (> 100 colonies = high level excretion, < 100 colonies = low excretion, or no colonies = no excretion) [3]. This cohort also provided an opportunity to follow the dynamics of the fecal microbiota of veal calves under real-life conditions by the collection of additional fecal samples. Characterization of the ESBL-producing *E. coli* will be published elsewhere.

Sampling began upon the arrival of batches of new calves in October and November 2015. A batch was defined in the study as a group of calves entering the farm at the same time and reared together until slaughter (Additional file 10: Fig. S7). Seven days after arrival, 15 calves were randomly selected from 50 on each farm and included in the study: five with high levels of ESBL-producing *E. coli* excretion, five with low-level excretion, and five with no excretion. The 15 calves selected per farm were then sampled by rectal swabbing bimonthly

for ESBL-producing *E. coli* excretion follow-up until departure of the batch to the slaughterhouse. To study the dynamics of the fecal microbiota, additional samples were collected at days 7 and 21 on the farms, then monthly for 5 months, for a total of seven samples per calf (Fig. 1). Swabs were placed immediately in portable coolers with ice packs, shipped to the ANSES lab in Lyon, France, and stored at -80 °C.

The calves were 14 days old when they arrived at the farms and were mainly fed milk replacer, which is reconstituted from cow milk powder with hot water, throughout fattening. Their diet was also supplemented with a small amount of solid feed from the first weeks. The daily quantity of milk powder was divided by the mean weight of a male Holstein calf of the corresponding age each day (assuming that the calves arrived at the farms at 14 days of age), which enabled us to estimate a proxy for the dose of milk powder consumed on these farms throughout the fattening period. Thus, the integrators' recommended doses of milk powder, of which approximately 45% consists of lactose, were estimated for farms B and C.

Collective antibiotic treatments were recorded by the three farms throughout fattening (Fig. 1). Antibiotics were always used at therapeutic doses and administered orally in water or milk replacer. All treatments upon entry to the fattening farms were set up to prevent gastrointestinal disorders, whereas treatments during the course of fattening were used to treat respiratory diseases. 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). On farm A, calves received a 10-day course of colistin and sulfonamides (first day to day 10), followed by another 10-day course of tetracycline (day 11 to day 20). Later during fattening, they received a one-day treatment of doxycycline on day 53 and a one-day treatment of amoxicillin on day 135. On farm B, calves received a one-day treatment of doxycycline and erythromycin on day 26, a one-day treatment of tetracycline on day 90, and a one-day treatment of amoxicillin on day 101. On farm C, calves received a six-day course of sulfonamides and trimethoprim (from day 3 to day 8) and a seven-day course of tetracycline (day 10 to day 16). They also received a five-day course of doxycycline (day 20 to day 24) and four-day courses of spiramycin (day 25 to day 28) and tetracycline (day 80 to day 83).

### DNA extraction, 16S rRNA gene sequencing, and Escherichia-specific quantitative PCR

Genomic DNA was extracted from rectal swabs using the DNEasy PowerSoil kit (QIAGEN, Venlo, Netherlands). The cotton tips of frozen swabs were broken off directly into bead tubes. The tubes were
incubated at 70 °C for 10 min, as previously described [48]. The remaining steps were performed according to the manufacturer's instructions, with an additional overnight incubation step with elution buffer at 4 °C. Extracted DNA was stored at -20 °C. The V4 region of the 16S rRNA gene from each sample was amplified using the primers 515fB (GTGYCAGCMGCCGCGGTAA) and 806rB (GGACTACNVGGGTWTCTAAT), modified to contain a barcode sequence between the primer and Illumina adaptor sequences, as previously described [49, 50]. Dual-barcoded libraries were sequenced on an Illumina MiSeq machine (MiSeq Reagent Kit V3, 600 cycles) according to the manufacturer's specifications to generate paired-end reads of 300 bases in length.

E. coli populations were quantified by qPCR targeting the 16S rRNA gene sequence specific to the *Escherichia* genus. The *Escherichia* genus is a good proxy of *E. coli* in cattle, as *Escherichia* cryptic clades and other species of the genus represent < 4% in non-human mammal feces [51]. We verified this assumption by plating three randomly selected swabs from distinct calves of our study on Drigalski plates. We determined the species of 50 colonies per plate using MALDI-TOF and Clermont typing PCR [52]. All colonies were confirmed to be *E. coli*. The sequences of the forward (CATGCCGCGT GTATGAAGAA) and reverse (CGGGTAACGT CAATGAGCAAA) primers and probe (FAM-TATTAA CTTTACTCCCTTCCTCCCCGCTGAA-TAMRA) were obtained from [53]. Each DNA sample (approximately 20 ng) was added to 30 µl PCR mixture containing 15 µl TaqMan Universal PCR master mix II 2X (Applied Biosystems, Life Technologies, Carlsbad, California, USA), 300 nM of each primer, 100 nM fluorescent probe, and bovine serum albumin at a final concentration of 0.1 µg/ µl (New England BioLabs, Evry, France), as previously described [54]. A standard curve was generated using known amounts of DNA of the archetypal ED1a *E. coli* strain for each experiment. Products were detected using an Applied Biosystems Prism 7500 instrument.

## Processing of 16S rRNA gene sequences

The quality of the reads was verified using FASTQC [55] and they were processed using mothur (version 1.35.1) [56, 57]. Contigs were generated by assembling forward and reverse reads. Low-quality contigs were discarded if the total length was outside 289 to 292 bases, if there were more than five ambiguous bases ("N"), or if homopolymer runs exceeded five bases. After a clean-up step, sequences were aligned to those of the SILVA reference database (February 2017, release 128) [58]. OTU assignment was made after clustering the sequences with a similarity cutoff of 97%. Singletons, duplicates, and triplicates were discarded. The taxonomy of each detected OTU was obtained using the RDP quality-

controlled, aligned, and annotated Bacterial and Archaeal 16S rRNA gene sequence database [59]. Chimeric sequences were removed after de novo chimera detection using the VSEARCH tool, version 2.3.4 [60]. Sequences flagged as chloroplasts, mitochondria, or eukaryotes were discarded from the dataset. We assessed the taxonomic composition of the microbiota, focusing on the phyla and genera, and focused on the five most abundant genera in each sample. For each calf, the OTUs detected in a sample were classified according to their detection in the previous sample to determine the monthly degree of change of the calf microbiota at the OTU level. For each time point, we determined the proportion of OTUs that were not previously detected and the proportion of OTUs that were not detected in the previous sample. This analysis was performed to identify concurrent patterns of acquisition and persistence of OTUs among calves.

A rarefaction step was performed before computation of the α- and β-diversity metrics. Several candidate rarefaction thresholds were assessed for several characteristics, such as the number of samples below the threshold, the number of samples for which the sampling effort would result in more than 25% and 50% of the total number of sequences, and the proportion of sequences sampled in the most abundant sequence sample. The threshold was set to 47,000 sequences, which was the best compromise between these four characteristics. Then, α-diversity and β-diversity metrics were computed from the rarefied samples. α-diversity is defined as the ecological diversity within samples, whereas β-diversity is defined as the dissimilarity between samples. The Shannon diversity index and the number of observed OTUs were computed as α-diversity metrics. Unweighted and raw-weighted Unifrac distances were computed as β-diversity metrics [61]. For each calf, distances were computed between consecutive samples, and the distances between all calves were computed at the first, second, and last sampling.

In summary, the available data for each sample were (i) the relative abundance of taxa at the phylum and genus level, (ii) the α-diversity metrics Shannon index and the number of observed OTUs, (iii) the weighted and unweighted Unifrac distances to the previous and next sample for the same calf, and (iv) the absolute number of *E. coli* per gram of feces. In addition, Unifrac distances were computed at the first, second, and last sampling between calves.

## Statistical analyses

### Comparing community structures between and within calves over time

We performed PERMANOVA tests [62] using weighted and unweighted Unifrac distance matrices to evaluate
the effects of time and the farm on the calf's microbiota. Tests were performed using 1000 permutations. We constrained the permutations within each calf to account for repeated measures. The weighted and unweighted Unifrac distances between calves at the first, second, and last samplings were represented in heatmaps. Moreover, the Unifrac distances between consecutive samples for each calf were represented by spaghetti plots to assess the temporal variability of its microbiota composition.

### Temporal modelling of α-diversity indices

We built linear mixed-effects models to study the temporal dynamics of microbiota diversity. One-slope and two-slope models were tested. For the two-slope model, the break was set at day 21 for the farms after visual inspection of the raw data. The equation of the two-slope model is shown below:

$$y_{i,j,k} = \theta_{0i,k} + \theta_{1i,k} \times t_{1i,j,k} + \theta_{2i,k} \times t_{2i,j,k} + \varepsilon_{i,j,k}$$

where $y_{i,j,k}$ is the observed Shannon index or the number of observed OTUs at the $j^{th}$ day of calf $i$ from farm $k$. $\theta_{0i,k}$, $\theta_{1i,k}$, $\theta_{2i,k}$ are the intercept, first slope, and second slope, respectively. $t_{1i,j,k}$ and $t_{2i,j,k}$ represent the time before or on the 21st day and after, respectively, and $\varepsilon_{i,j,k}$ the residual error.

The farm effect was introduced for each parameter and the calves set as random effects. The LRT was used to compare the fit between the candidate models. For each parameter of the final model, farms were grouped when the effects were not significantly different. We assumed that the random effects and residual errors were independent and had a normal distribution, with a mean of 0. Evaluation of the final model was conducted using basic goodness-of-fit plots.

### Temporal modelling of the absolute number of E. coli per gram of feces

First, we performed Spearman's correlation between the number of bacteria from the *Escherichia* genus, estimated by quantitative PCR, and the relative abundance of the *Escherichia* genus to assess the consistency of the two techniques. Second, we built linear mixed-effects models to study the temporal dynamics of the absolute number of *E. coli* and tested polynomial functions of time. The farm effect was introduced for each parameter and the calves set as random effects. The equation of the quartic model is shown below:

$$y_{i,j,k} = \theta_{0i,j,k} + \theta_{1i,j,k} \times t_{1i,j,k} + \theta_{2i,j,k} \times t_{2i,j,k} + \theta_{3i,j,k} \times t_{3i,j,k} + \theta_{4i,j,k} \times t_{4i,j,k} + \varepsilon_{i,j,k}$$

where $y_{i,j,k}$ is the observed number of *E. coli/g* and $\theta_{0i,j,k}$, $\theta_{1i,j,k}$, $\theta_{2i,j,k}$, $\theta_{3i,j,k}$, $\theta_{4i,j,k}$ the coefficients of each term of the polynomial function of time; $t_{1i,j,k}$ the $j^{th}$ day of calf $i$, and $\varepsilon_{i,j,k}$ the residual error.

The number of random effects was reduced by a backward approach using the Akaike Information Criterion (AIC), starting with the random effect of the coefficient with the highest degree. Selection of the final model and grouping by the farm effect were performed as for the α-diversity indices. We also made the same assumptions of normality and independence and evaluated the final model using basic goodness-of-fit plots.

### Determination of the influence of antibiotic treatment on temporal predictions of α-diversity indices and the absolute number of E. coli

Alterations of fecal microbiota, such as selective depletion of bacterial populations and the reduction of ecological diversity, are generally regarded as gut microbiota dysbiosis markers following antibiotic treatment [63–65]. We investigated such effects on the temporal dynamics of the calves' fecal microbiota by considering that dysbiosis linked to antibiotic treatment could have occurred for all samples for which a collective antibiotic treatment was given within 15 days before sampling and antibiotic treatment still ongoing or not at the time of sampling (Fig. 1).

The influence of antibiotic treatment on the intercept of α-diversity indices and the absolute number of *E. coli* models was tested by introducing a covariable representing the existence of treatment within the 15 days before sampling. Samples were collected in such a time window on two dates at farm A (days 7 and 21), one date at farm B (day 106), and four dates at farm C (days 7, 21, 35, and 91, Fig. 1).

### Exploring associations between the abundance of genera and the dose of milk powder

A link between the relative abundance of genera and the dose of milk powder was explored for farms B and C. The Spearman correlation test was used to look for positive associations between the dose of milk powder and the relative abundance of genera. As multiple tests were performed, the P values were adjusted using the Bonferroni correction. We searched for the presence of the β-galactosidase enzyme, which hydrolyzes lactose to the monosaccharides glucose and galactose, in the available annotated genomes in the NCBI genome and protein databases for genera that were found to have a significant moderate to strong positive association with the dose of milk powder ($r_s > 0.4$) [37]. We also blasted the protein sequence of the LacZ β-galactosidase enzyme in the NCBI protein database using the program blastp [36, 37].

As *E. coli* is a lactose-fermenting species, we searched for an association between the estimated daily doses of milk powder and the predicted number of *E. coli/g* for farms B and C using Spearman's correlation test. The
farm-predicted numbers of *E. coli/g* were obtained by adding the estimates of the farm parameters to the intercept in our final model.

Means were presented with standard deviations. All statistical analyses were performed using R software (R version 3.1.0) [66]. Mixed effect models were built using nlme [67] and PERMANOVA tests were performed using the “adonis” function in the vegan package [68].

**Supplementary information**

Supplementary information accompanies this paper at https://doi.org/10.1186/s42523-020-00052-6.

**Additional file 1 Fig. S1.** Sequence and OTU distributions after bioinformatics processing. (a) Distribution of the number of 16S rRNA gene V4 region sequences in samples after quality filtering. (b) Distribution of the number of OTUs in samples after clustering sequences with a similarity cutoff of 97%. The inner lines in the boxplots represent the median, the edges show the first and third quartiles, and the whiskers extend to the 5th and 95th percentiles in (a) and (b). (c) Rarefaction curves for 16S rRNA gene V4 region sequences. Each curve corresponds to a sample. The red vertical line represents the chosen rarefaction threshold.

**Additional file 2 Fig. S2.** Heatmaps of the β-diversity unweighted Uni-frac distances matrix for the (a) first sampling (day 7), (b) second sam-pling (day 21), and (c) last sampling (day 161 for farms A and B and day 147 for farm C). Yellow squares indicate low Unifrac distances, whereas dark red squares indicate high Unifrac distances. Calves are ordered ac-cording to farms in both lines and columns. The means ± standard devia-tions for each sampling on each farm are shown in the lower triangles.

**Additional file 3 Fig. S3.** Observed intra-calf β-diversity weighted Uni-frac distances between consecutive samplings for (a) farm A, (b) farm B, and (c) farm C. The dots indicate the Unifrac distances between consecu-tive samples from the same calf.

**Additional file 4 Fig. S4.** Relative abundance of the five most abundant taxa at the genus level for all calves throughout the fattening period. For each panel, the first and second days represent the sampling date for farm C and farms A and B, respectively. Relative abundance of the five most abundant taxa are given for (a) days 35 and 49, (b) days 63 and 77, (c) days 91 and 106, and (d) days 119 and 133. Other detected taxa are depicted by the white bars. Calf IDs are provided at the top of the panels and are ordered according to farms. The color scale of the dots beneath the bar graphs represents the distribution of the Shannon index values. The color key refers to the phylum of each taxa and each palette was built to maximize the distinctiveness between shades.

**Additional file 5 Table S1.** OTUs detected in or absent from previous samples and shared by calves over time. The sampling, ranges of proportions of calves, and OTU taxonomies are presented in each layer. The lists of OTUs consist of (a) OTUs simultaneously detected in two consecutive samples in more than 25% of the calves, (b) OTUs that were not previously detected and that simultaneously appeared in more than 25% of the calves, and (c) OTUs that were simultaneously lost by more than 25% of the calves.

**Additional file 6 Tables S2, S3, and S4.** Tables S2, S3, and S4 contain the estimated parameters for the final models of the Shannon index, the number of observed OTUs, and the absolute number of *E. coli/g*, respectively.

**Additional file 7 Fig. S5.** Dynamics of the mean observed and predicted number of observed OTUs for each farm. Predicted dynamics of the number of observed OTUs, without and with the antibiotic-treatment effect, in the final model are represented in panels (a) and (b), respectively. The mean values ± standard deviations of the observed data for each farm are represented by the dashed bars. Model-predicted profiles and their 95% confidence bands are represented by the solid lines and bands, respectively. Antibiotic treatments during sampling or within

15 days before sampling are color-coded by farm and indicated above the x-axis in panel (b).

**Additional file 8 Table S5.** List of the Spearman's rank correlation coefficients between genera found in calves' feces and the dose of milk powder estimated for farms B and C. The genera with a significant positive correlation and a Bonferroni corrected *p*-value < 0.05 are highlighted in green.

**Additional file 9 Fig. S6.** Relative abundance of the genera *Megasphaera*, *Enterococcus*, *Dialister*, and *Mitsuokella* as a function of the dose of milk powder. Each point represents a sample. These four genera had the highest significant positive correlation with the estimated dose of milk powder in farms B and C. Values on the x-axis correspond to samples in which the corresponding genus was not detected by 16S rRNA gene sequencing.

**Additional file 10 Fig. S7.** Veal calves on fattening farms (a) on the first day, corresponding to 14 days of age, and (b) at 115 days of age, during the third month of fattening.

**Additional file 11 Table S6.** Absolute number of *Escherichia coli* per gram of feces in all samples, estimated by *Escherichia*-specific quantitative PCR.

**Abbreviations**

AIC: Akaike Information Criterion; ANSES: French Agency for Food, Environmental and Occupational Health & Safety; ESBL: Extended spectrum β-lactamase; LRT: Likelihood ratio test; OTU: Operational taxonomic unit; PERMANOVA: Permutational multivariate analysis of variance; qPCR: Quantitative PCR; rRNA: Ribosomal ribonucleic acid

**Acknowledgements**

We are grateful to Antoine Bridier-Nahmias, Hervé Perdry, and David Gordon for helpful discussions and Mélanie Magnan, Caroline Wybraniec, and Gilles Collin for technical assistance.

**Authors' contributions**

MH, JYM, FM, and ED conceived and designed the study. MH and JYM collected the samples. MM performed the laboratory assays and carried out the bioinformatics analyses. MM, TTN, and FM carried out the statistical analyses of the data. MM generated the figures. MM and ED wrote the manuscript. MH, TTN, JYM, and FM revised and edited the draft. The authors read and approved the final manuscript.

**Funding**

This work was supported in part by grants from La Fondation pour la Recherche Médicale to M.M. (grant number FDM20150633309) and E.D. (équipe FRM 2016, grant number DEQ20161136698), and from the European Union's Horizon 2020 Research and Innovation Programme under Grant Agreement No. 773830 (Project ARDIG, EJP One Health) and INTERBEV (Protocol N° SECU-15-31). The funding organizations were not involved in the study design, or collection, analysis, or interpretation of the data, or writing of the manuscript.

**Availability of data and materials**

Sequencing reads were deposited as entire raw data in the European Nucleotide Archive repository (ENA) under the BioProject ID PRJEB33072, separately for each sample. The code used for processing the 16S rRNA gene sequences was originally developed by Kozich et al. [57] and is available at the website https://www.mothur.org/wiki/MiSeq_SOP (access on January 2017). The unrarefied OTU table and the corresponding taxonomic classification analyzed during the current study are available from the corresponding author upon reasonable request. The *Escherichia*-specific qPCR dataset has been included as Additional file 11: Table S6.

**Ethics approval and consent to participate**

This study was declared to the CNIL, the French office responsible for protecting personal data, supporting innovation, and preserving individual liberties. No further ethical approval was needed since this study did not involve any experimentation on animals (only rectal swabs were sampled) and since we did not collect and register any personnel opinion of the participants.
**Consent for publication**

Not applicable.

**Competing interests**

The authors declare that they have no competing interests.

**Author details**

¹Université de Paris, IAME, INSERM, Site Xavier Bichat, 16 rue Henri Huchard, F-75018 Paris, France. ²Unité Antibiorésistance et Virulence Bactériennes, Université de Lyon - ANSES, Laboratoire de Lyon, Lyon, France. ³AP-HP, Hôpital Bichat-Claude Bernard, Département d'Épidémiologie, Biostatistiques et Recherche Clinique, F-75018 Paris, France. ⁴AP-HP, Hôpital Bichat-Claude Bernard, Laboratoire de Génétique Moléculaire, F-75018 Paris, France.

Received: 12 April 2020 Accepted: 4 September 2020

Published online: 15 September 2020

**References**

1. Sans P, de Fontguyon G. Veal calf industry economics. Rev Médecine Vét. 2009;160:420-4.

2. Jarrige N, Cazeau G, Morignat E, Chanteperdrix M, Gay E. Quantitative and qualitative analysis of antimicrobial usage in white veal calves in France. Prev Vet Med. 2017;144:158-66.

3. Gay E, Bour M, Cazeau G, Jarrige N, Martineau C, Madec J-Y, et al. Antimicrobial usages and antimicrobial resistance in commensal *Escherichia coli* from veal calves in France: evolution during the fattening process. Front Microbiol. 2019;10:792.

4. Pardon B, Catry B, Dewulf J, Persoons D, Hostens M, De Bleecker K, et al. Prospective study on quantitative and qualitative antimicrobial and anti-inflammatory drug use in white veal calves. J Antimicrob Chemother. 2012; 67:1027-38.

5. Ji S, Jiang T, Yan H, Guo C, Liu J, Su H, et al. Ecological restoration of antibiotic-disturbed gastrointestinal microbiota in foregut and hindgut of cows. Front Cell Infect Microbiol. 2018;8:79.

6. Holman DB, Yang W, Alexander TW. Antibiotic treatment in feedlot cattle: a longitudinal study of the effect of oxytetracycline and tulathromycin on the fecal and nasopharyngeal microbiota. Microbiome. 2019;7:86.

7. Oultram J, Phipps E, Teixeira AGV, Foditsch C, Bicalho ML, Machado VS, et al. Effects of antibiotics (oxytetracycline, florfenicol or tulathromycin) on neonatal calves' faecal microbial diversity. Vet Rec. 2015;177:598.

8. Haley BJ, Kim S-W, Salaheen S, Hovingh E, Van Kessel JAS. Differences in the microbial community and resistome structures of feces from preweaned calves and lactating dairy cows in commercial dairy herds. Foodborne Pathog Dis. 2020;17:494-503.

9. Smith HW, Crabb WE. The faecal bacterial flora of animals and man: its development in the young. J Pathol Bacteriol. 1961;82:53-66.

10. Meale SJ, Li S, Azevedo P, Derakhshani H, Plaizier JC, Khafipour E, et al. Development of ruminal and fecal microbiomes are affected by weaning but not weaning strategy in dairy calves. Front Microbiol. 2016;7:582.

11. Kolenda R, Burdukiewicz M, Schierack P. A systematic review and meta-analysis of the epidemiology of pathogenic *Escherichia coli* of calves and the role of calves as reservoirs for human pathogenic *E. coli*. Front Cell Infect Microbiol. 2015;5:23.

12. Litvak Y, Byndloss MX, Tsolis RM, Bäumler AJ. Dysbiotic *Proteobacteria* expansion: a microbial signature of epithelial dysfunction. Curr Opin Microbiol. 2017;39:1-6.

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

14. 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.

15. Wickramasinghe HKJP, Anast S-ES, Serão NVL, Appuhamy JADR.N. Beginning to offer drinking water at birth increases the species richness and the abundance of *Faecalibacterium* and *Bifidobacterium* in the gut of preweaned dairy calves. J Dairy Sci. 2020;103(5):4262-74.

16. 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.

17. Klein-Jöbstl D, Schornsteiner E, Mann E, Wagner M, Drillich M, Schmitz-Esser S. Pyrosequencing reveals diverse fecal microbiota in Simmental calves during early development. Front Microbiol. 2014;5:622.

18. Malmuthuge N, Griebel PJ, Guan LL. The gut microbiome and its potential role in the development and function of newborn calf gastrointestinal tract. Front Vet Sci. 2015;2:36.

19. Yeoman CJ, Ishaq SL, Bichi E, Olivo SK, Lowe J, Aldridge BM. Biogeographical differences in the influence of maternal microbial sources on the early successional development of the bovine neonatal gastrointestinal tract. Sci Rep. 2018;8:1-14.

20. Badman J, Daly K, Kelly J, Moran AW, Cameron J, Watson I, et al. The effect of milk replacer composition on the intestinal microbiota of pre-ruminant dairy calves. Front Vet Sci. 2019;6:371.

21. Jami E, Israel A, Kotser A, Mizrahi I. Exploring the bovine rumen bacterial community from birth to adulthood. ISME J. 2013;7:1069-79.

22. Dill-McFarland KA, Breaker JD, Suen G. Microbial succession in the gastrointestinal tract of dairy cows from 2 weeks to first lactation. Sci Rep. 2017;7:40864.

23. Malmuthuge N, Li M, Fries P, Griebel PJ, Guan LL. Regional and age dependent changes in gene expression of toll-like receptors and key antimicrobial defence molecules throughout the gastrointestinal tract of dairy calves. Vet Immunol Immunopathol. 2012;146:18-26.

24. Foditsch C, Pereira RVV, Siler JD, Altier C, Warnick LD. Effects of treatment with enrofloxacin or tulathromycin on fecal microbiota composition and genetic function of dairy calves. PLoS One. 2019;14:e0219635.

25. Xie G, Duff GC, Hall LW, Allen JD, Burrows CD, Bernal-Rigoli JC, et al. Alteration of digestive tract microbiome in neonatal Holstein bull calves by bacitracin methylene disalicylate treatment and scours. J Anim Sci. 2013;91: 4984-90.

26. Song Y, Malmuthuge N, Steele MA, Guan LL. Shift of hindgut microbiota and microbial short chain fatty acids profiles in dairy calves from birth to pre-weaning. FEMS Microbiol Ecol. 2018;94(3):10.1093/femsec/sec179.

27. 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.

28. Liu J, Taft D, Gomez M, Johnson D, Treiber M, Lemay D, et al. The fecal resistome of dairy cattle is associated with diet during nursing. Nat Commun. 2019;10(1):4406.

29. Bosman AB, Wagenaar JA, Wagenaar J, Stegeman A, Vernooij H, Mevius D. Quantifying antimicrobial resistance at veal calf farms. PLoS One. 2012;7:e44831.

30. Pantophlet AJ, Gilbert MS, van den Borne JJGC, Gerrits WJ, Roelofsen H, Priebe MG, et al. Lactose in milk replacer can partly be replaced by glucose, fructose, or glycerol without affecting insulin sensitivity in veal calves. J Dairy Sci. 2016;99:3072-80.

31. Sabarly V, Bouvet O, Glodt J, Clermont O, Skumik D, Diancourt L, et al. The decoupling between genetic structure and metabolic phenotypes in *Escherichia coli* leads to continuous phenotypic diversity. J Evol Biol. 2011;24:1559-71.

32. Ozawa A, Freter R. Ecological mechanism controlling growth of *Escherichia coli* in continuous flow cultures and in the mouse intestine. J Infect Dis. 1964;114:235-42.

33. Fabich AJ, Leatham MP, Grissom JE, Wiley G, Lai H, Najar F, et al. Genotype and phenotypes of an intestine-adapted *Escherichia coli* K-12 mutant selected by animal passage for superior colonization. Infect Immun. 2011;79: 2430-9.

34. Sabarly V, Aubron C, Glodt J, Balliau T, Langella O, Chevret D, et al. Interactions between genotype and environment drive the metabolic phenotype within *Escherichia coli* isolates. Environ Microbiol. 2016;18:100-17.

35. Wu S, Zhang Y, Lu R, Xia Y, Zhou D, Petrof EO, et al. Intestinal epithelial vitamin D receptor deletion leads to defective autophagy in colitis. Gut. 2015;64:1082-94.

36. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403-10.

37. NCBI Resource Coordinators. Database resources of the National Center for Biotechnology Information. Nucleic Acids Res. 2018;46(D1):D8-D13.

38. Parks D, Chuvochina M, Waite D, Rinke C, Skarshewski A, Chaumeil P-A, et al. A standardized bacterial taxonomy based on genome phylogeny substantially revises the tree of life. Nat Biotechnol. 2018;36(10):996-1004.

39. Beukers AG, Zaheer R, Goji N, Amoako KK, Chaves AV, Ward MP, et al. Comparative genomics of *Enterococcus spp* isolated from bovine feces. BMC Microbiol. 2017;17:52.
40. Fortina MG, Ricci G, Mora D, Manachini PL. Molecular analysis of artisanal Italian cheeses reveals *Enterococcus italicus* sp. nov. Int J Syst Evol Microbiol. 2004;54:1717-21.

41. Svec P, Vancanneyt M, Koort J, Naser SM, Hoste B, Vihavainen E, et al. *Enterococcus devriesei* sp. nov., associated with animal sources. Int J Syst Evol Microbiol. 2005;55:2479-84.

42. Collins MD, Jones D, Farrow JAE, Kilpper-Balz R, Schleifer KH. *Enterococcus avium* nom. Rev., comb. nov.; *E. casseliflavus* nom. Rev., comb. nov.; *E. durans* nom. Rev., comb. nov.; *E. gallinarum* comb. nov.; and *E. malodoratus* sp. nov. Int J Syst Evol Microbiol. 1984;34:220-3.

43. Trosvik P, Stenseth NC, Rudi K. Convergent temporal dynamics of the human infant gut microbiota. ISME J. 2010;4:151-8.

44. Smith HW. The development of the bacterial flora of the faeces of animals and man: the changes that occur during ageing. J Appl Bacteriol. 1961;24: 235-41.

45. Alipour MJ, Jalanka J, Pessa-Morikawa T, Kokkonen T, Satokari R, Hynönen U, et al. The composition of the perinatal intestinal microbiota in cattle. Sci Rep. 2018;8:10437.

46. Vandeputte D, Kathagen G, D'hoe K, Vieira-Silva S, Valles-Colomer M, Sabino J, et al. Quantitative microbiome profiling links gut community variation to microbial load. Nature. 2017;551(7681):507-11.

47. Hordijk J, Wagenaar JA, van de Giessen A, Dierixk C, van Essen-Zandbergen A, Veldman K, et al. Increasing prevalence and diversity of ESBL/AmpC-type β-lactamase genes in *Escherichia coli* isolated from veal calves from 1997 to 2010. J Antimicrob Chemother. 2013;68:1970-3.

48. Costello EK, Lauber CL, Hamady M, Fierer N, Gordon JI, Knight R. Bacterial community variation in human body habitats across space and time. Science. 2009;326(5960):1694-7.

49. Apprill A, McNally S, Parsons R, Weber L. Minor revision to V4 region SSU rRNA 806R gene primer greatly increases detection of SAR11 bacterioplankton. Aquat Microb Ecol. 2015;75:129-37.

50. Parada AE, Needham DM, Fuhrman JA. Every base matters: assessing small subunit rRNA primers for marine microbiomes with mock communities, time series and global field samples. Environ Microbiol. 2016;18:1403-14.

51. Smati M, Clermont O, Bleibtreu A, Fourreau F, David A, Daubié A-S, et al. Quantitative analysis of commensal *Escherichia coli* populations reveals host-specific enterotypes at the intra-species level. MicrobiologyOpen. 2015;4: 604-15.

52. Clermont O, Christenson JK, Denamur E, Gordon DM. The Clermont *Escherichia coli* phylo-typing method revisited: improvement of specificity and detection of new phylo-groups. Environ Microbiol Rep. 2013;5:58-65.

53. Penders J, Thijs C, Vink C, Stelma FF, Snijders B, Kummeling I, et al. Factors influencing the composition of the intestinal microbiota in early infancy. Pediatrics. 2006;118:511-21.

54. Smati M, Clermont O, Le Gal F, Schichmanoff O, Jauréguy F, Eddi A, et al. Real-time PCR for quantitative analysis of human commensal *Escherichia coli* populations reveals a high frequency of subdominant phylogroups. Appl Environ Microbiol. 2013;79:5005-12.

55. Andrews S. FastQC: a quality control tool for high throughput sequence data. 2010. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc.

56. Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol. 2009;75:7537-41.

57. 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:5112-20.

58. Pruesse E, Quast C, Knittel K, Fuchs BM, Ludwig W, Peplies J, et al. SILVA: a comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB. Nucleic Acids Res. 2007;35:7188-96.

59. Cole JR, Wang Q, Cardenas E, Fish J, Chai B, Farris RJ, et al. The ribosomal database project: improved alignments and new tools for rRNA analysis. Nucleic Acids Res. 2009;37:D141-5.

60. Rognes T, Flouri T, Nichols B, Quince C, Mahé F. VSEARCH: a versatile open source tool for metagenomics. PeerJ. 2016;4:e2584.

61. Lozupone C, Knight R. UniFrac: a new phylogenetic method for comparing microbial communities. Appl Environ Microbiol. 2005;71:8228-35.

62. Anderson MJ. A new method for non-parametric multivariate analysis of variance. Austral Ecol. 2001;26:32-46.

63. Palleja A, Mikkelsen KH, Forslund SK, Kashani A, Allin KH, Nielsen T, et al. Recovery of gut microbiota of healthy adults following antibiotic exposure. Nat Microbiol. 2018;3:1255-65.

64. Korpela K, Salonen A, Virta LJ, Kekkonen RA, Forslund K, Bork P, et al. Intestinal microbiome is related to lifetime antibiotic use in Finnish pre-school children. Nat Commun. 2016;7:10410.

65. Burdet C, Nguyen TT, Saint-Lu N, Sayah-Jeanne S, Hugon P, Sablier-Gallis F, et al. Change in bacterial diversity of fecal microbiota drives mortality in a hamster model of antibiotic-induced Clostridium difficile colitis. Antimicrob Agents Chemother. 2018;63:382.

66. R Development Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for statistical computing; 2005.

67. Pinheiro J, Bates D, DebRoy S, Sarkar D, R Core Team. Nlme: linear and nonlinear mixed effects models. R package version 3.1-137; 2018.

68. Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, et al. Vegan: community ecology package. 2019.

## Publisher's Note

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

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

<figure><img src="image_10.png" /></figure>