Animal housing conditions were described in a previously published paper (Monger et al., 2022). The experimental design was approved by the Animal Protection Committee of Université Laval prior to the experiments (2019057–1). Briefly, six Yorkshire-Landrace male pigs at 55 days of age were transported to the Université Laval animal research facility. At 66 days, they underwent surgery to implant an ileal-T canula, as previously described (Wubben et al., 2001). They were allowed a 11-day recovery period before the beginning of the treatments, as it was demonstrated that inflammation at the cannulation site was negligible 7 days after the surgery (Zhang et al., 2014). The six pigs received three different treatments: macrolide antibiotic (tylvalosin TYL; 250 g/ton of feed; Aivlosin^® (17% tylvalosin), ECO Animal Health Princeton, NJ, USA), probiotic Pediococcus acidilactici MA18/5M (PA; 10⁸ CFUs/day; Biopower^® PA, Lallemand Animal Nutrition, Milwaukee, WI, USA), and a combination of both antibiotic and probiotic treatment at the same concentrations as the individual treatments. Each pig received the three treatments in a different order, and a period of recovery was given between each treatment with a control diet not supplemented, according to a crossover design (Figure 1). The periods of treatment and recovery were 3 weeks each.

Digesta and faeces from each pig were collected at the end of each treatment period for a total of 36 faeces samples and 35 digesta samples. The digesta from pig number two at the Ctl2 period could not be collected as there was no digesta at the time of sampling. All the samples were stabilized using the PERFORMAbiome•GUT | PB-200 sampling kit according to the manufacturer’s instructions (DNAgenotek, Ottawa, Ontario, Canada) and as suggested elsewhere (Monger et al., 2022).

The DNA was extracted from the samples with the QIAamp PowerFecal Pro DNA Kit (QIAGEN, Toronto, Canada) using 250 µL of stabilized faeces or digesta. The DNA was quantified using the PicoGreen kit (Invitrogen, Waltham, Massachusetts, USA). Integrity of the DNA was assessed by electrophoresis on an agarose gel.

Short chain fatty acids (SCFAs) were quantified as previously described (Laforge et al., 2023). Briefly, samples were stored at −80°C, then thawed and homogenized in water. The homogenized samples were centrifugated and SCFAs were extracted by liquid-liquid extraction. The SCFAs were subjected to gas chromatography coupled with a flame ionization detector (GC-FID Shimadzu, Kyoto, Japan) following the protocol published by Roussel et al. (2022).

Shotgun sequencing libraries were generated from 50 ng of gDNA using the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England BioLabs, Whitby, Canada) as per the manufacturer’s recommendations. Sequencing was done on an Illumina NovaSeq 6000 by the Genome Quebec Centre of Expertise and Services (Montréal, Canada). Base calling was performed with Illumina’s Real Time Analysis software version 3.4.4. A minimum of 70 million paired reads was generated for each of the samples. The optical duplicates in the reads were removed with the clumpify tool from bbmap version 38.96 (Bushnell, 2014). The reads were then filtered using fastp version 0.23.1 (Chen S. et al., 2018), and the sequences were mapped on pig genome sequences using bowtie2 version 2.4.4 (Langmead and Salzberg, 2012) and samtools version 1.17 (Li et al., 2009) to keep only those of the microbiome. The reference genome used was Sscrofa11.1 (RefSeq: GCF_000003025.6). The number of reads at each step is presented in Supplementary Table S1. The whole genome shotgun sequencing data set was deposited in the NCBI Sequence Read Archive database under the BioProject ID PRJNA1049315.

The reads were then analyzed and coassembled with the SqueezeMeta pipeline, version 1.5.1, using default parameters (Tamames and Puente-Sanchez, 2019). Only contigs with a length greater than or equal to 200 bp were retained for further analysis. The SQMtools R package (Puente-Sanchez et al., 2020) was used to import the data in R. The vegan R package (Dixon, 2003) was used for alpha diversity analysis, and the phyloseq (McMurdie and Holmes, 2013) and microbial (Kai, 2021) packages were used for beta diversity and linear discriminant analysis (LDA) analysis, respectively. The cutoffs for LDA were a score greater than or equal to two and an adjusted p-value less than 0.05. The pvclust package was used for hierarchical clustering of the data with 10,000 bootstrap replications (Suzuki and Shimodaira, 2006).

The assembly was screened for plasmidic contigs using PlasForest version 1.4 (Pradier et al., 2021). The plasmid sequences were reconstructed from contigs using the MOB-suite version 3.1.4 (Robertson and Nash, 2018). ARGs were identified by mapping the reads against the Comprehensive Antibiotic Resistance Database (CARD) (Alcock et al., 2020) using MetaProtMiner with a minimum similarity of 90% and a query cover of 80% (Galiot et al., 2023).

The 16S rRNA gene was sequenced and analyzed as described elsewhere (Laforge et al., 2023). Briefly, the v3–v4 region of the 16S DNA was sequenced at the Plateforme d’analyse génomique of the Institut de Biologie Intégrative et des Systèmes (Université Laval, Quebec City, Canada). The sequences were pretreated and taxonomically assigned with the R package DADA2 (Callahan et al., 2016). During filtration, the first 17 nucleotides of the forward reads and the first 21 nucleotides of the reverse reads were trimmed to remove primers. Dereplication, sample inference, chimera identification, and merging of the paired-end reads were performed using the default parameters. Taxonomic assignment was conducted using the SILVA rRNA database (release 138.1) with the naive Bayesian classifier method (utilizing the assignTaxonomy command of the DADA2 package). The sequencing data set was deposited in the NCBI Sequence Read Archive database under the BioProject ID PRJNA1049315.

Statistical analyses were performed using the prism software version 9.2.0. For multiple comparisons, data normality was assessed using the Shapiro-Wilk test. A Kruskal–Wallis test followed by a post hoc Dunn’s test was used to analyze alpha diversity and short-chain fatty acids. A Tukey test was used to compare antibiotic resistance gene scores. Adjusted p-values were employed to evaluate the significance of the results. Statistical analysis of the beta diversity was performed with the adonis function from the vegan package, which is a function to compute permanova tests with 10,000 iterations. Correlation analysis of the SCFAs data (lactic, propionic, butyric, and acetic acids) with bacterial taxa was done with R using the ccrepe package and 10,000 iterations (Schwager et al., 2019). Only correlations with compositionality corrected p-values less than 0.05 were considered significant.

An antibiotic treatment (tylvalosin), a probiotic treatment (P. acidilactici), and a combination of both treatments were administered to six pigs to investigate the resilience and adaptation of their gut microbiome using metagenomic analysis of their faeces and ileal digesta samples. To assess the effect of the treatments on the microbial diversity within the samples (alpha diversity), the Shannon and Simpson diversity indexes were computed for all samples and the treatments were compared (Figure 2). Both indexes are a measure of diversity considering evenness and number of taxa; the Shannon index is more influenced by richness, and the Simpson index is more influenced by evenness of repartition and less influenced by rare species.

No significant difference was observed (p > 0.05) between treatments for either digesta (Figures 2A, C) or faeces (Figures 2B, D) samples. The metataxonomic analysis of 16S rRNA gene sequences did not reveal any significant difference between treatments for either sample type (Supplementary Figure S1).

The dissimilarity between the samples (beta diversity) was visualized by a principal coordinate analysis (PCoA) plot of the Bray-Curtis distance (Figure 3). The treatments had no significant effect on the microbiota composition of the digesta (p > 0.05; Figure 3A). In the faeces (Figure 3B), all groups were different from the first control period (all p < 0.0096). No other groups were found to be significantly different. As expected, the digesta and the faeces were significantly different from each other (p < 0.0001, Figure 3C) and the variability between faeces samples was greater than the variability between digesta samples. Also, 16S rRNA gene metagenomics validated that the microbiota remained stable between treatments (Supplementary Figure S2).

The microbiota composition is similar between treatments (Supplementary Figure S3) and as expected (Holman et al., 2017). The digesta is predominantly composed of the Lactobacillaceae and Peptostreptococcaceae families, accounting for more than half of the microbial composition. Notably, there is a significant increase in Streptococcaceae observed in the samples taken after probiotic treatment. The composition of the feces is more evenly distributed among various taxa, including several families of unclassified Eubacteriales, Bacteroidales, Firmicutes, and Clostridia, as well as the Prevotellaceae family. However, as anticipated, the fecal microbiome is initially largely colonized by bacteria from the Prevotellaceae family at the beginning of the experiment, followed by a subsequent reduction.

Marker taxa were identified using LDA and some of the bacteria varied significantly (p < 0.05) between the antibiotic and probiotic groups (Supplementary Figure S4). Interestingly, the genera that were significantly different between these two groups were only detected during treatment with the probiotic. The most amplified genus in the digesta was the Methanosphaera archaea, while in the faeces it was an unclassified bacterial genus belonging to the Akkermansiaceae family.

For each treatment period, the ARG score (calculated by MetaProtMiner) of the previous recuperation period was subtracted from the score calculated for each treatment period to assess only the effect of the treatment (Figures 4A, B). There was no significant difference in digesta samples between treatments (p > 0.05, Figure 4A). However, in the faecal samples, the ARG score change caused by the treatments was significantly higher for the antibiotics and combined treatments than for the probiotic treatment (Figure 4B).

The ARG score for the preceding treatment was subtracted from the score for each recuperation period to verify the resilience of the animals’ microbiome following treatment. During the recuperation periods following treatments, the ARG score change was not significantly different between treatments for the digesta (Figure 4C). However, in the faeces, ARGs decreased during the recuperation period following the antibiotic and combined treatments (Figure 4D). The change during recuperation from the combined treatment was significantly lower than during recuperation from the probiotic treatment (p = 0.04). There was a tendency (p < 0.061) for the decrease in ARGs to be smaller during recuperation after antibiotic treatment than during recuperation from probiotic treatment. In the faeces for the combined treatment, the decrease in ARGs during recuperation was smaller than the increase during treatment (p = 0.02), while the difference was not significant for the antibiotic treatment (p = 0.31).

Sample grouping based on the resistome profile was done with a PCoA plot using Bray-Curtis dissimilarity (Figure 5). In the digesta (Figure 5A), significant differences were observed between the initial microbiome (control) and the antibiotics treatment (p = 0.046), but the comparison was not significant between the initial microbiome and the combined treatment, although the p-value is close to 0.05 (p = 0.051). Differences were also found between the probiotic treatment and antibiotic treatment (p = 0.037) and the combined treatment (p = 0.005). In the faeces (Figure 5B), all groups were different (all p < 0.003) except for the combined and antibiotic treatments (p = 0.7534). The resistome profile was also different between the total samples coming from the digesta and the faeces (p < 0.0001, Figure 5C). More variability in the resistome profile was observed in the faeces samples from the probiotic treatment compared with the combined and antibiotic treatments. Interestingly, the opposite trend was observed in the digesta samples: the samples from the probiotic treatment group were less variable than the combined and antibiotic treatment groups. Globally, there was also more variability in the faeces samples than in the digesta samples (Figure 5C).

Marker resistance genes between the probiotic and antibiotic treatments were identified with LDA analysis (Supplementary Figure S5); no markers were found between the combined and antibiotic treatments. A heatmap of the abundance score of those genes across samples is shown in Figure 6. In both faeces and digesta, hierarchical clustering analysis of the transitions in the whole resistome between the control and treatments showed that the changes induced by the probiotic treatment clustered together, while the changes induced by the antibiotic or combined treatments clustered together (Figure 6). The ARGs that were considered markers with LDA were amplified in the transition to antibiotic and combined treatments and were stable or diminished in the transition to probiotic treatment. Five ARG markers were found in the digesta (Figure 6A), of which ermB, ermG, and ermQ were macrolide resistance genes. A total of 14 ARG markers were found for the faeces samples (Figure 6B), and five of these (ermB, ermG, ermT, ermQ, and ermF) were macrolide resistance genes. This is not surprising, as tylvalosin is an antibiotic of the macrolide class. Interestingly, all the genes identified as markers in the digesta samples were also markers in the faeces samples.

After screening of the total assembly for plasmidic sequences, 38,565 of 9,444,776 contigs were identified to be putatively plasmidic. These sequences were compared with the CARD database, and 27 contigs were carrying at least one ARG and two contigs were carrying two ARGs. No contigs were found carrying more than two ARGs. Plasmid sequences were reconstructed from contigs, and those carrying ARGs are described in Table 1. Two plasmids, AA998 and AA448, were carrying multiple ARGs, and both were considered as conjugative plasmids by the tool mob_suite. The plasmid AA448 is associated with Klebsiella pneumonia and was carrying genes for resistance against the tetracycline, gentamycin, and sulfonamide classes of antibiotics. The plasmid AA998 was associated with Salmonella enterica serovar Typhi and was carrying resistance genes against the beta-lactam, tetracycline, chloramphenicol, and macrolide classes of antibiotics (Table 1).

To identify marker plasmidic contigs between the antibiotic and probiotic treatments, LDA analysis was performed (Supplementary Figure S6). In the digesta, all maker plasmids were amplified in the probiotic treatment; however, none of the marker contigs were carrying ARGs (Supplementary Figure S6). In the faeces, all markers for plasmidic contigs were amplified in the antibiotic treatment. Four of the contigs that were markers in the faeces were carrying ARGs, and three of these belonged to plasmids that were either conjugative or mobilizable.

Samples were grouped based on the transitions in the whole plasmidome between the control and treatments using hierarchical clustering with heatmaps of the transition of the ARG-carrying plasmids (Supplementary Figure S7). The whole plasmidome of the samples did not cluster according to treatments.

To more specifically assess the role of the treatments on functional profile, the relative abundance of functions, as defined by the Pfam database (Finn et al., 2008), was established and markers were identified using LDA analysis. In the digesta, all the marker functions were more frequent in the probiotics group and no horizontal transfer functions were identified as markers (Figures 7A). In the faeces, two of the functions identified as markers for antibiotic treatments were recombinase (PF07508 [Recombinase]) and resolvase (PF00239 [Resolvase, N terminal domain]); these two functions play a role in gene exchange (Figures 7B) (Johnson and Grossman, 2015; Rice, 2015).

Four SCFAs were measured in the samples: acetic acid, lactic acid, propionic acid, and butyric acid; these are important metabolites for host health (Bergman, 1990; De Vadder et al., 2014; Jacobson et al., 2018; Schulthess et al., 2019). The SCFA content was compared between treatments, and there were no significant differences in the quantity of SCFAs measured between treatments in the digesta (Supplementary Figure S8) or the faeces (Supplementary Figure S9).

The individual SCFA contents were also compared between the digesta and the faeces (Figure 8). The contents of acetic, propionic, and butyric acid were higher in the faeces (p < 0.0001) while lactic acid was higher in the digesta (p < 0.0001).

The correlations between the SCFA content and bacterial species were also tested in the faeces and the digesta. In the digesta, various bacteria, including many Clostridium species, were negatively correlated with lactic acid and positively correlated with the other SCFAs. In contrast, lactobacilli were positively correlated with lactic acid, and negatively or not correlated with other SCFAs (Figure 9A). In the faeces, lactic acid did not seem to follow any trend, but there was a clear separation between bacteria that were positively correlated with acetic acid and negatively with butyric and propionic acids, and bacteria that were positively correlated with butyric and propionic acid and negatively correlated with acetic acid (Figure 9B).