The animal trial was reviewed and approved by the ethical committee of the Faculty of Veterinary Medicine and Bioscience Engineering of Ghent University (EC 2020-016). The care and treatment of the birds were in full compliance with the national and European legislation concerning animal welfare (European Commission., 2010; Royal Decree., 2013). Throughout the study, feed and drinking water were provided ad libitum.

Eighty-four healthy 1-day-old male broiler chicks (Ross 308), with similar body weights (bws) of approximately 55 g were randomly allocated to one of four treatments [n = 7 cages (3 chickens each) per treatment]: (1) standard commercial feed (Table 1, Farm mash 1&2, Versele-Laga, Deinze, Belgium) for 37 days (control group), 2) antibiotic-treated group 1 [OPT, 12.5 mg/kg bw/day Baytril^® 10% oral solution (Bayer, Diegem, Belgium) consecutively via drinking water from day 21 to day 23 of age], (3) antibiotic-treated group 2 (HIGH, 50 mg/kg bw/day Baytril^® 10% oral solution consecutively via drinking water from day 21 to day 23 of age), and (4) synbiotic supplemented group (OPT_PS), in which chickens received the commercial feed supplemented with a synbiotic product (PoultryStar^®, BIOMIN GmbH, Austria) during their lifetime, in conjunction with receiving enrofloxacin at a dose of 12.5 mg/kg bw/day via drinking water from day 21 to day 23 of age. The age of the birds was selected to reflect the real commercial setting and scenario for antimicrobial use, since colibacillosis frequently occurs in the second half of the broiler growing period. Additionally, the 3-day treatment duration was in compliance with the label of Baytril. The synbiotic supplement was composed of 2 × 10⁸ CFU/g multi-species probiotic: Enterococcus faecium (1.20 × 10⁸ CFU/g), Bifidobacterium animalis (6 × 10⁷ CFU/g), Lactobacillus salivarius (2 × 10⁷ CFU/g), and prebiotic inulin. The amount of mL enrofloxacin that needed to be added in the drinking water was calculated using the following formula (Equation 1) (Temmerman et al., 2021):

Immediately before (day 20) and after enrofloxacin treatment (day 24), and 2 weeks after antibiotic withdrawal (day 37), one bird per replicate was sacrificed and digesta samples from the cecum were collected in Eppendorf cups, immediately snap frozen (liquid N2), and stored at −20°C. The animals were euthanized with pentobarbital IV (sodium pentobarbital 20%, Kela, Hoogstraten, Belgium) dosed at approximately 100 mg/kg.

DNA was extracted using the QIAamp PowerFecal Pro DNA Kit (Qiagen, Antwerp, Belgium), which as shown in our recent studies, has minimum effects on the microbial community structure and promising results in terms of the DNA integrity for chicken cecum samples (Wegl et al., 2021). The purity of the DNA was assessed by evaluating the 260/280 and 260/230 absorbances using NanoDrop^® (ND-1000, Thermo Fisher Scientific, Merelbeke, Belgium), while Quantus^® (Promega, Leiden, The Netherlands) was used to estimate the DNA concentration. Finally, DNA integrity was visually validated by gel electrophoresis. The DNA was sent to LGC Genomics GmbH, Berlin, Germany, where paired-end sequencing libraries were constructed using the Illumina Nextera XT Library Preparation Kit (Illumina Inc., San Diego, CA) followed by sequencing on the Illumina NextSeq 500 platform using high-output chemistry (2 × 150 bp) according to the manufacturer’s instructions.

After demultiplexing and adaptor trimming (trimmomatic v. 0.39) (Bolger et al., 2014) taxonomic profiles of the demultiplexed reads and taxa relative abundance estimates were generated using Kraken2 (Wood et al., 2019) and Bracken (Lu et al., 2017), with kraken2 confidence set to 0.05 and bracken threshold to t −10. The abundances of ARGs were quantified by mapping the reads against the hand-curated AMR database MEGARes v2.0 (Doster et al., 2020) using USEARCH (v10) (Edgar, 2010). In MEGARes 2.0, the nodes of the acyclic hierarchical ontology included four antimicrobial compound types, 57 classes, 220 mechanisms of resistance, and 1,345 gene groups that classified the 7,868 accessions. High confidence matches to the sequence in the MEGARes database were obtained by considering the entire coverage of the query reads against ARG genes with an identity threshold of 90% (parameters were set as “-usearch-global -id 0.9, maxaccepts 1, threads 50”) (Ghanbari et al., 2019). The variance explained by each variable included in the trial was quantified by omnibus testing on Bray–Curtis dissimilarity matrices from metagenome data, using PERMANOVA with the adonis function in the R package Vegan (Dixon, 2003). Ordination was performed with log-transformed normalized reads on two dimensions with the “prcomp” ordinate function from the package “stats” v3.6.2 and visualized using the “fviz_pca” function from the “factoextra” function from the FaceFacefactoextra package v1.0.7.

Linear regression analysis was performed considering random differences between the treatment groups at the first sampling (for the control and enrofloxacin treatments) as covariate terms in the model if the richness (the number of unique taxa) and inverse Simpson diversity (the number and relative abundance of unique taxa) were affected by the different treatments. The Inverse Simpson diversity index was selected due to its biological interpretation aspect and since it does not tend to be as affected by sampling effort as other diversity indexes. The Mann–Whitney test was used for statistical comparisons between the pairs of groups at each sampling time point. Log-fold changes in abundance (of taxa and ARGs) between groups were determined by a negative binomial generalized linear model using DESeq2 version 1.26.0. (Love et al., 2014) in R, considering the statistical significance at BH-corrected P ≤ 0.10. In all statistical analyses, the cage was considered as the experimental unit.

The average weights per treatment group (7 animals per group, 28 in total) and associated variability at the end of the animal trial are shown in Figure 1. The broiler weight ranges at day 37 were 1.25–1.95, 1.37–2.20, 1.25–2.00, and 1.58–2.00 kg for the control, HIGH, OPT, and OPT_PS groups, respectively. The corresponding averages are 1.59, 1.82, 1.54, and 1.83 kg, respectively. Chickens that received the highest enrofloxacin dose (HIGH group) and the optimized enrofloxacin dose supplemented with the synbiotic (OPT_PS group) had higher average weights at 37 days of age when compared to the control group, although only the chickens in the OPT with synbiotic supplementation (OPT_PS) group showed a much higher homogeneity in weight gain and a significant difference in weight group (p = 0.034, Figure 1). Interestingly, the difference between the optimized dosing (OPT) and OPT_PS was statistically significant (p = 0.026), suggesting an important role of synbiotic application in restoring physiological balance after perturbation by antibiotics and, consequently, improving the weight gain.

Shotgun metagenomics of the 84 cecum samples generated over 600 million paired-end sequence reads with a mean number of eight million reads per sample (see Supplementary File 1, which contains detailed information on the samples, including library size). On average, 23% of the reads per sample were classified to the species level using Kraken2 with a confidence threshold of 5%. Using the MEGARes v2.0 database with a 90% gene cut-off fraction, 1,281,069 reads were aligned to 1,000 antimicrobial resistance genes across treatments. Permutational Multivariate Analysis of Variance (PERMANOVA) statistical analysis calculated on Bray–Curtis dissimilarity matrices from metagenome data showed that the age of the birds (sampling time points: days 20, 24, and 37) was the largest driver of variation in both cecum microbiota and ARG resistome profiles, with antibiotic and synbiotic application capturing less variation (Figure 2A). In line with this, principal component analysis (PCA) revealed a clear separation of the sample based on the day of sampling (age) (Figures 2B,C). Regarding the effect of antibiotic administration on global diversity, the lower, optimized dosage of enrofloxacin (OPT) resulted in a higher variation compared to the higher dosage (HIGH) on the composition and structure of the taxonomical (Adonis R2, 5.9% vs. 3.6%) and the ARG resistome profiles (Adonis R2, 3.9% vs. 3.1%) in broiler chickens. Synbiotic application (OPT_PS) significantly explained 3.1 and 4.6% of the variation (Adonis R2) of the cecum microbiota and ARG resistome profiles, respectively (Figure 2A).

Further, pairwise post hoc comparisons¹ showed a significant difference in the profiles of the relative abundance of taxa and ARGs (with the exception of the optimized dosage effect on the taxa profiles) between the enrofloxacin (both HIGH and OPT) and control treatment at day 24 of sampling. No significant differences were seen for the sample taken at day 20 or for those taken 2 weeks after the withdrawal of enrofloxacin at day 37 (Table 2), indicating the resilience of the bacterial communities carrying ARGs from antibiotic perturbation. On the other hand, OPT_PS showed a significant difference with the control groups in the relative abundance of taxonomy (with the exception of day 20) and resistome profiles at days 20 and 37 (Table 2). On day 24, immediately after an optimized dosage of enrofloxacin was applied in combination with the synbiotic product, the microbiota and ARG resistome profiles of the samples clearly diverged from the control group (P = 0,006), suggesting a synergistic effect on the chicken cecal metagenome (Table 2).

Linear regression analysis, taking into account the random differences between the high and optimized doses of enrofloxacin with control treatments on day 20, revealed a significant association of enrofloxacin application with the reduction of the microbiota Observed index from days 20 to day 24 (HIGH vs. control coefficient estimate day 24 = −41.143 with P = 2.55e^–05, OPT vs. control coefficient estimate day 24 = −33.857 with P = 0.000414, and OPT_PS vs. control observed index coefficient estimate day 24 = −41.857 with P = 1.91e^–05). Two weeks after the withdrawal of the antibiotic, although the microbial community richness recovered to some extent (HIGH vs. control CE day 37 = −26.429 with P = 0.00507, OPT vs. control EC day 37 = −26.000 with P = 0.005795, and OPT_PS vs. control CE day 37 = −20.143, with P = 0.030770), the cecum microbiota Observed index was still significantly lower in all groups than in the control. For OPT_PS at day 20, the observed index for the cecum microbiota was higher compared to those in the control (187 ± 9 vs. 183 ± 6), although the difference was not statistically significant.

Calculation of the Inverse Simpson diversity index revealed that while the cecum microbiota diversity in the animals in the control group increased (not significantly) from day 20 (inverse Simpson 11.248 ± 1.819) to day 24 of age (12.842 ± 2.617), for the broiler chickens receiving the HIGH and OPT dosage, there was either no change (HIGH group, IS 12.048 ± 4.256–12.210 ± 4.0189) or a reduction (not significant; OPT group, IS 11.909 ± 3.059–10.536 ± 3.024) in microbial diversity (Figure 3A). For OPT_PS, the reduction in diversity was less pronounced (Figure 3A). At day 20 of age, OPT_PS showed generally higher diversity compared to the control (IS 12.960 ± 3.259 vs. 11.248 ± 1.819, respectively). While the withdrawal of enrofloxacin application did not seem to have a large impact on the recovery of the diversity for the animals in the HIGH and OPT_PS groups, it resulted in a significant increase in the cecum microbiota diversity in the optimized dosage group (Figure 3A).

A significant association of enrofloxacin application with the change in the diversity of ARGs was revealed [linear regression analysis, F_{(11, 72)} = 2.258, p = 0.02016 with an R2 of 0.26]. Over the course of the enrofloxacin application, that is, days 21–23, the diversity of ARGs in the cecum metagenome of the birds receiving enrofloxacin significantly increased (P < 0.05), regardless of dosage (Figure 3B). While OPT_PS had the lowest ARG diversity at day 20, the application of the optimized dosage of enrofloxacin on the top of the synbiotic supplement resulted in the largest increase (P < 0.01) on the diversity of ARGs from days 20 to 24 compared to all other treatments (Figure 3B). Yet, the cecum ARG diversity of birds in the OPT group (17.393 ± 0.727) at day 20 was significantly higher than that in the control (16.477 ± 1.651, P = 0.620), HIGH (16.176 ± 1.412, P = 0.097), and OPT_PS (16.089 ± 0.617, P = 0.011) groups, suggesting an alleviation effect of synbiotic application on the expansion of broiler cecal resistome (Figure 3B). Withdrawal of antibiotic application resulted in a pronounced reduction in ARG diversity, with OPT_PS showing significantly lower diversity of ARGs compared to the control group (P = 0.011) at day 37 (Figure 3B).

Taxonomic profiling of the cecum samples in terms of the most abundant microbial phyla and species over the course of the study is shown in Figure 4. In total, 1,573 different microbial species were observed, of which species belonging to Firmicutes, followed by species belonging to Actinobacteria and Proteobacteria, were the most dominant members of the cecal microbiota (Figure 4).

Differential abundance (DA) analysis of the species-level taxonomic assignments using the DESeq2 tool revealed a significant effect of enrofloxacin application on shaping the cecal microbiota over the course of antibiotic application. Among the different groups of the impacted taxa, species belonging to the Proteobacteria phylum were among the most affected taxa due to enrofloxacin application followed by species from Firmicutes and Actinobacteria (Figure 5A and Supplementary Data 1). In fact, over the course of the enrofloxacin application (days 21–23), both high and optimized dose application resulted in a significant (FDR-corrected p-value < 10%) increase in the abundance of 17 proteobacterial species, while significantly reducing the abundance of 15 proteobacterial species (Figure 5B and Supplementary Data 1). Two weeks after the last day of antibiotic application, the number of affected microbial species belonging to the Proteobacteria decreased, with only four and five species showing a significantly higher abundance in the HIGH and OPT groups compared to the control, respectively (Figure 5B and Supplementary Data 1). Synbiotic application significantly (FDR-corrected p-value < 10%) reduced the abundance of 15 proteobacterial species, including Escherichia coli and Campylobacter spp., on day 20 (Supplementary Data 1). Interestingly, the OPT_PS group, showed a lower number of significantly increased Proteobacteria species compared to the OPT group (11 vs. 17, Figure 5B and Supplementary Data 1). Furthermore, with regard to the Firmicutes phylum, while species from the Lactobacillus, Bacillus, and Anaerobutyricum showed a significant increase (FDR-corrected p-value < 10%) in the cecum microbiota of chickens in OPT_PS compared to the control, several microbial pathogenic species such as Enterococcus cecorum, Enterococcus hirae, and Enterococcus gallinarum were significantly decreased (FDR-corrected p-value < 10%) (Supplementary Data 1). At day 37, 2 weeks after the withdrawal of antibiotics, the highest number of significantly decreased Proteobacteria was observed in the broiler chickens in OPT_PS (Figure 3).

Antibiotic resistome profiling of the cecal metagenome showed that, regardless of antibiotic medication, all samples harbored a diverse range of ARGs (Figure 6A). Genes encoding MLS (mean = 29.64 ± 0.045 across all samples) and tetracycline (Mean = 27.70 ± 0,051 across all samples) resistance were the most predominant classes of ARGs to which the reads were aligned, followed by oxazolidinone and aminoglycoside resistance genes (Figure 6A). Correlation analysis showed that among the different dominant groups of the microbial phyla, Proteobacteria had the highest positive correlation (Spearman correlation, R = 0.7, P = 3.1e^–13, Figure 6B) with the diversity of ARGs in the cecal metagenome of the broiler chickens (Figure 7B), with Actinobacteria showing a weak positive correlation (Spearman correlation R = 0.29, P = 0.0085) and Firmicutes showing a strong negative correlation (Spearman correlation R = −0.6, P = 1.2e^–10).

In line with the observation that several members of Proteobacteria significantly increased in the cecal metagenome of broiler chickens receiving enrofloxacin, DA analysis revealed a significant increase in the abundance of ARGs encoding resistance to antibiotics belonging to different classes of antibiotics (mainly MLS, aminoglycosides, and tetracyclines) in these two groups from days 20 to 24, with the optimized dosage application resulting in a twofold higher number of DA ARGs compared to high dose application (32 DA ARGs vs. 17 DA ARGs, respectively; Figures 7, 8 and Supplementary Data 1). Examples of the increased ARGs in the samples from the HIGH and OPT treatments taken at day 24 include cat and catA (encoding chloramphenicol acetyltransferases), ermB (encoding 23S rRNA methyl-transferases), and MLS23S (macrolide-resistant 23S rRNA mutation) (Figures 7B,C). The most strongly diminished ARGs were the asmA (multi-drug RND efflux regulator) and ptsL (fosfomycin target mutation) genes, with an approximately −6 log2 fold change (Figures 7B,C). Other decreased genes included those encoding gyrA, pbp2 (Penicillin binding protein2), and CTX (class A beta-lactamases), as well as tetA, tetBP, and tet32 (Figures 7B,C). Withdrawal of antibiotic administration seems to result in the recovery of ARG resistome to the same level as the control group, as there were almost no DA ARGs in the cecum metagenomes of broiler chickens between the high/optimized dosage and those in the control group (Figure 8). For the samples taken from the OPT_PS group at day 20, where only synbiotic products were applied, the abundance of 13 ARGs was significantly lower than that in the control group, with three ARGs showing a significantly higher abundance in the chicken cecum metagenomes (Figures 7A, 8 and Supplementary Data 1). At day 24, for OPT_PS, an alleviating effect of synbiotic application on the ARG burden was observed (Figures 7D, 8 and Supplementary Data 1). Fifty-three ARGs (vs. 32 ARGs in OPT) were significantly lower in abundance in this group than in the control group, and only 13 genes (vs. 32 ARGs in OPT), mainly belonging to tetracycline resistance, showed significantly higher abundance compared to the control group (Figure 7D and Supplementary Data 1). No DA ARGs were found in the animals in this group at day 37 compared to the control group.