In total, samples from 30 mice passed quality controls, 15 in each group. One animal from the donor B group at day 0 (D0) was excluded due to low microarray signal intensity values. Four hundred fifty-four OTUs (operational taxonomic units—see the “Methods” section for the definition as we used the dada2 pipeline [22]), were identified following filtering and were assigned to 51 genera (mean sequences per sample following filtering 6089, range 4219–8835, rarefaction curves presented in Additional file 1: Figure S1A).

Baseline microbiome composition modulates susceptibility and recovery following antibiotics exposure

Donor A recipient mice (donor A mice) had a markedly different community compared to mice from donor B. Donor A mice had a microbiota composition dominated by Prevotella (mean 26.6% (SD 13.6%)) and Faecalibacterium (mean 24.9% (SD 12.5%)), with Bacteroides (mean 17.3% (SD 7.8%)) being the third most dominant genus (Additional file 2). In donor B mice, the dominant genera were Bacteroides (mean 37.9% (SD 12.2%)) and Parabacteroides (mean 19.6% (SD 7.6%)), with less than 1% mean abundance of both Faecalibacterium and Prevotella (Fig. 1b, c). Samples for each donor were aggregated at each time point and submitted to enterotype clustering with the MetaHIT data set [23, 24] (following manual filtering to assure all genera matched) (Additional file 1: Figure S1B). This confirmed clustering of donor A mice with the Prevotella enterotype and donor B mice with the Bacteroides enterotype. Online reference-based enterotyping [25] similarly classified these samples to the Prevotella and Bacteroides enterotypes.

The response to antibiotic administration was markedly different between the two donor groups. Immediately following co-amoxiclav treatment (D8—day 8), donor A mice demonstrated an increase in Prevotella sequences from baseline (D0) (from mean 27.4% (SD 5.9%) to mean 44.9% (SD 3.7%)) with a significant reduction in this genus by late recovery (mean 13.9% (SD 7.7%)) and a corresponding increase in Faecalibacterium OTUs at these time points (D0—mean 13.1% (SD 10.5%), D8—mean 21.3% (SD 5%), D11 (day 11)—mean 23.6% (SD 13%), D18 (day 18)—mean 38.5% (SD 7.4%); Fig. 1b, Additional file 2). There was a trend of decreasing proportional abundance in Bacteroides throughout the study in the donor B group, although this did not reach significance at any time point (D0—49.3% (SD 17.5%), D8—40.9% (SD 9.7%), D11—31.9% (SD 12%), D18—32.2% (SD 4%)), while Parabacteroides remained largely stable throughout (Additional file 2). Diversity was not significantly altered by antibiotic treatment (Fig. 1d, Additional file 3). Comparing different time points using DESeq2 for donor A mice, there was a significant reduction in Clostridium XIVb at D8 from baseline (D0), a significant reduction in Roseburia at D8 from D0, a significant reduction in Prevotella between D8 and D18, a significant increase in Blautia at D18 compared to D8, and a significant increase in Dorea at D18 compared to D8 (Fig. 2a).

Fig. 2 a Genera differentially abundant (adjusted P value < 0.01) in donor A mice by DESeq2. b Similar to a but for donor B mice. c Weighted unifrac distance with per group. Results of the corresponding two-way PERMANOVA are presented in Table 1. d PCoA of weighted unifrac distances for donor A alone, with genera that correlate significantly with PCoA axes (Spearman’s correlation, P value after correction 0.1) and PERMANOVA R2 and P value for univariate comparison. e PCoA for donor B as in d with weighted unifrac distances. No significant correlations were present. PERMANOVA R2 and P value for univariate comparison Full size image

In contrast to donor A mice, only one genus, Ruminococcus, increased significantly in donor B mice at D18 (Fig. 2b). Overall, there was a much more marked disruption of the community structure in mice from donor A mice when compared to donor B mice, evident in PCoA plots of weighted unifrac distance (Fig. 2c). To determine if there was a significant difference in terms of the effect of donor and day post antibiotics, we performed a two-way PERMANOVA, which indicated a significant interaction between donor and day (Table 1). A post hoc analysis demonstrated differences only for donor A, at pairwise contrasts between D0 vs D18, D8 vs D11, and D8 vs D18 (false discovery rate (FDR) P values < 0.05). Correlation of the PCoA scores with genera demonstrated a correlation between a number of key taxa, including Prevotella and Faecalibacterium, associated with the changes in donor A mice (Fig. 2d). No correlations were evident for donor B mice (Fig. 2e). While we did not have a no-antibiotic control group, we designed the study such that antibiotic-induced changes to the microbiota must occur independently in two cages. To ensure this was the case, we used a distance-based method to determine if within-cage differences over time (Additional file 4: Figure S2A) were significantly different from time point-specific distances, including between-cage effects at D8 (post antibiotics) (Additional file 4: Figure S2B and S2C). We observed a significant shrinking in the distance post antibiotics in donor A mice but did not observe a similar reduction in donor B mice, suggesting that the effects of antibiotics in this group were not above the baseline within-cage variability over the course of the study (Additional file 4: Figure S2C).

Table 1 Results of two-way PERMANOVA of weighted unifrac distances (random seed for rooting tree, 1234). P values (***< 0.001, *< 0.05) Full size table

Taken together, these results demonstrated a marked effect of co-amoxiclav administration on the microbiota of donor A mice, with notable shifts in the dominant genus Prevotella and a number of less abundant genera. In contrast, the microbiota profile of donor B mice remained relatively unaffected.

Antibiotics induce dramatic changes on host transcriptome, which overlap independently of baseline microbiota composition

Gene expression was assessed initially with PCA for all mice combined (Fig. 3a). This demonstrated a pattern associated with antibiotic treatment, with the main separation between early (D0 and D8) and recovery (D11 and D18) time points. Figure 3b represents differentially expressed genes between D18 and D0 (1845 in total), which again demonstrate clear clustering between early (D0 and D8) and recovery (D11 and D18) time points. Gene Ontology (GO) analysis with the goana function in limma identified a large number of significant GO terms between different contrasts, so the top 100 were submitted to REVIGO and of the resulting output, those with a − log10(P value) > 5 were plotted. This analysis demonstrated significant changes in many terms, prominently biological processes concerned with extracellular matrix organization, cellular adhesion, signaling, cell migration, developmental processes, and circadian rhythm (Fig. 3c).

Fig. 3 a PCA plots of gene expression data following antibiotic treatments for all mice combined, clustered by day with respect to antibiotic treatment. b Heatmap of differentially expressed genes at the D18 vs D0 contrast. c GO enrichment analysis, followed by submission to REVIGO, stratified by contrast Full size image

A combined analysis of both groups in our study also demonstrated significant differential gene expression alterations at D8 relating to genes associated with circadian rhythm, consistent with previously reported findings following antibiotic treatment, including altered expression of Per genes, Cipc, Prkaa2, NR1D2/Rev-ErbA, Ciart, NPAS2/CLOCK, NFIL3/E4BP4, and ARNTL/BMAL1 [26] (Additional file 5: Figure S3). Recent work has precisely explored the relationship between circadian rhythm and the microbiome and identified diurnal biogeographic changes in microbial proximity to the epithelium which orchestrate epithelial and host circadian processes [27,28,29].

When the same analysis was repeated for each donor group individually, more significantly differentially expressed genes were detected in the donor A than donor B group (Fig. 4a–c). While only a limited number of differentially expressed genes were evident for donor B mice, there was notable overlap in those genes that were with donor A mice (Fig. 4c). Additionally, heatmaps of contrasts D11 vs D8 (Fig. 4d) and D18 vs D8 (Fig. 4e) including these genes, while dominated by donor A, still demonstrated clear clustering based on early (D0, D8) and late (D11, D18) time points, again suggesting that even though these genes did not meet the threshold for differential expression for donor B, there was a common pattern in response to co-amoxiclav. GO terms with at least four genes were submitted to REVIGO (donor A, Fig. 4f; donor B, Fig. 4g) demonstrated common activation of circadian rhythm terms at D8 compared to recovery time points (D11 and D18), while terms relating to immune system processes, cell adhesion and response to stimulus were evident in donor A mice at D18 when compared to the early time points (D0 and/or D8).

Fig. 4 Gene expression analysis for both donor groups individually. a PCA as in Fig. 3a, isolating donor A gene expression. b PCA for donor B expression alone. c Contrasts with differential genes and the relative contribution from donor A and donor B. Red represents donor A, blue donor B, and white overlap. d D11 vs D8 heatmap of differential gene expression. e D18 vs D8 heatmap of differential gene expression. f GO enrichment analysis following submission to REVIGO, stratified by contrast for donor A. g As in f but for donor B Full size image

Taken together, these results demonstrate a marked effect of co-amoxiclav administration on the colonic transcriptome of humanized mice. Gene expression analysis identified changes in genes involved in extracellular matrix organization, cellular adhesion and motility, developmental processes, circadian rhythm, and immune system processes. While the pattern of gene expression was similar regardless of the baseline (D0) microbiota, the degree of change appeared to be more marked in mice humanized by donor A than by donor B. As these were the mice that had a more marked alteration in the microbiota composition in response to co-amoxiclav, we next investigated whether correlation between the transcriptome and microbiota was evident.

Host transcriptome-microbiota correlation reveals marked covariance between the microbiota and GO pathways in donor A mice but not in donor B mice

Due to the large number of differentially expressed genes across the different contrasts, we used Gene Set Variation Analysis (GSVA) to identify gene pathways associated with antibiotic treatment. To perform correlation with microbial taxa, we used Hierarchical All-against-All significance testing [30], a recently developed technique for determining relationships between multi-omics datasets (Fig. 5, Additional files 6 and 7). This method identifies significant correlations between individual features of both datasets, as well as identifying correlations between clusters of features. Performed using Spearman’s correlation values on the genus-level abundance from each donor and the corresponding pathway abundance from GSVA, a large number of correlations were evident for donor A mice following multiple-hypothesis correction, including a number of dominant and/or differentially abundant bacterial taxa from that donor (Prevotella, Dorea, Faecalibacterium, Blautia, and Roseburia). Figure 5 represents a network plot of significant correlations between individual pathways and bacterial genera for donor A mice. The pathways are colored by high-level GO assignments. Furthermore, pathway-bacteria correlations, including additional clusters determined by HAllA, are represented in Additional file 6 (data) and Additional file 7: Figure S4. In contrast, no significant correlations between pathway expression and the microbiota were evident for donor B mice.

Fig. 5 Network plot of significant correlations between GO pathways drawn from GSVA and genera abundance for donor A mice determined by the HAllA method. Corresponding GO pathways are provided in Additional file 6 Full size image

Taken together, these results demonstrate correlation between host pathways and microbiome composition was unique to donor A mice. However, as the antibiotic-induced changes in donor B mice were considerably weaker, such correlations are harder to detect and higher numbers of mice would probably be required in this case.