Sow treatment induced long term changes in the nasal microbiota of the offspring
Globally, a total of 14,655,905 sequences were obtained by 16S amplicon (V3-V4) sequencing, with 12,084 different ASVs. Diversity analysis was calculated at a 12,124 read depth per sample, which required the elimination of one sample from the dataset due to low sequencing throughput (6,832 reads in a sample from the IP group at D7).
The effect induced by ceftiofur treatment in piglets’ nasal microbiota was evaluated when administered to sows before farrowing. Alpha diversity was estimated longitudinally (Chao and Shannon indexes), from birth to fattening at four time points. Initially (D0), alpha diversity was not different in piglets born to treated or non-treated sows, but piglets born to treated sows showed a higher inter-individual variability (Fig. 1). The treatment on sows induced a temporal increase in the alpha diversity of the nasal microbiota of piglets at D7 (Shannon diversity index P = 0.028) compared with the control group (Fig. 1). This rise in alpha diversity was not maintained through time and no significant differences were observed between TS and the control group at D21 or D49. When the beta diversity was analyzed, significant differences were detected at all sampling times, by either quantitative (Bray Curtis, weighted Unifrac) or qualitative (Jaccard, unweighted Unifrac) indexes. A stronger divergence was observed at D7 (Fig. 2), since the percentage of explanation measured on weighted Unifrac between groups was higher (R2 = 35%) at this timepoint compared to D21 (R2 = 28%) and D49 (R2 = 22%)
Ceftiofur treatment directly administered to piglets in early life (TP) produced similar trends in alpha diversity than the treatment applied to their sows (TS). A significant increase in alpha diversity (P = 0.009, Shannon diversity index) was detected at D7 in TP in comparison to the control group (Fig. 1). As observed in TP group, this increment was not maintained through time, and it was not observed at weaning (D21) or at the end of the study (D49). When beta diversity of TP and control groups was estimated, we detected differences at D7 and D21 using quantitative (Fig. 3) and qualitative metrics. At D49, the difference between these two groups was only detected using Jaccard distance (qualitative measurement). Importantly, the effect was more evident at D7 (Adonis function, R2 = 0.45, P = 0.012) than at the other timepoints (Adonis function, R2 = 0.24, and R2 = 0.15, at D21 and D49 respectively), underlying the different microbial composition in TP compared to the control soon after the antibiotic treatment.
Lastly, both treated groups (TP and TS) were compared to determine if the changes led by the antibiotic treatment were similar when administered to the piglets or to their sows. No significant differences were observed in alpha diversity at any sampling point. When beta diversity was analyzed, Jaccard, unweighted Unifrac and Bray Curtis metrics detected significative differences at D7 (P = 0.01, P = 0.026 and P = 0.03, respectively), but this difference was not observed using weighted Unifrac (Additional file 1: Fig. S1). Quantitative differences (weighted Unifrac and Bray Curtis) were not maintained at D21, while they were still detected with the qualitative analysis (Jaccard and unweighted Unifrac). Similar results that those found at D7 were obtained at D49 (unweighted Unifrac P = 0.006; Jaccard P = 0.008; Bray Curtis P = 0.012). Differences in the significance of quantitative and qualitative metrics pointed out that these differences were provoked by bacterial communities with low relative abundance.
To understand the microbial changes leading to these differences, we performed an analysis of the compositional microbiome with bias correction (ANCOM-BC) at each timepoint based on either amplicon sequence variants or at different taxonomic levels. At D7, we found 488 ASV in TS and 205 ASVs in TP that were differentially abundant when compared to the control group (ASVs; Additional file 2: Table S1), which was in agreement with the higher alpha diversity observed in this group. A significant increase in the relative abundance of ASVs from Bacteroidales (3.972% in TS, 3.534% in TP versus 0.357% in control) and Clostridiales (18.487% in TS, 7.733% in TP versus 2.452% in control) was observed at D7 in both antibiotic-treated groups. On the other hand, typical members of the nasal microbiota were decreased in the treated groups with respect to the non-treated control group, including Glaesserella, Moraxella and Rothia (Table 5).
Table 5
Relative abundance of genera belonging to the core microbiota from control group at weaning from each group at D7 and D21.
| D7 | D21 |
Taxa | Control | TS | TP | Control | TS | TP |
p__Proteobacteria;o__Pasteurellales;f__Pasteurellaceae;g__Glaesserella | 0.0183 | 0.0601 | 0.0017 | 0.2754 | 0.1094 | 0.1490 |
p__Bacteroidetes;o__Flavobacteriales;f__[Weeksellaceae];g__Bergeyella | 0.0000 | 0.0250 | 0.0000 | 0.1063 | 0.2005 | 0.0679 |
p__Proteobacteria;o__Pseudomonadales;f__Moraxellaceae;g__Enhydrobacter | 0.1252 | 0.0013 | 0.0010 | 0.0876 | 0.1265 | 0.1616 |
p__Firmicutes;o__Lactobacillales;f__Streptococcaceae;g__Streptococcus | 0.0243 | 0.0249 | 0.0145 | 0.0320 | 0.0746 | 0.0368 |
p__Proteobacteria;o__Caulobacterales;f__Caulobacteraceae;g__Caulobacter | 0.1005 | 0.0773 | 0.1593 | 0.0761 | 0.0801 | 0.1341 |
p__Proteobacteria;o__Pseudomonadales;f__Moraxellaceae;g__Moraxella | 0.3518 | 0.0500 | 0.0476 | 0.0673 | 0.0298 | 0.0190 |
p__Proteobacteria;o__Rhizobiales;f__Rhizobiaceae;g__Rhizobium | 0.0495 | 0.0307 | 0.0705 | 0.0402 | 0.0718 | 0.0785 |
p__Firmicutes;o__Lactobacillales;f__Lactobacillaceae;g__Lactobacillus | 0.0174 | 0.0530 | 0.1236 | 0.0022 | 0.0095 | 0.0052 |
p__Actinobacteria;o__Actinomycetales;f__Micrococcaceae;g__Rothia | 0.0069 | 0.0232 | 0.0100 | 0.0173 | 0.0026 | 0.0058 |
p__Bacteroidetes;o__Bacteroidales;f__Bacteroidaceae;g__Bacteroides | 0.0008 | 0.0087 | 0.0214 | 0.0025 | 0.0018 | 0.0044 |
p__Proteobacteria;o__Rhizobiales;__;__ | 0.0495 | 0.0307 | 0.0705 | 0.0003 | 0.0009 | 0.0017 |
p__Bacteroidetes;o__Bacteroidales;f__Prevotellaceae;g__Prevotella | 0.0008 | 0.0090 | 0.0064 | 0.0012 | 0.0015 | 0.0010 |
p__Firmicutes;o__Clostridiales;f__Ruminococcaceae;g__ | 0.0008 | 0.0046 | 0.0010 | 0.0020 | 0.0019 | 0.0009 |
p__Proteobacteria;o__Pseudomonadales;f__Moraxellaceae;g__Acinetobacter | 0.0151 | 0.0996 | 0.0816 | 0.0018 | 0.0048 | 0.0151 |
p__Proteobacteria;o__Rhizobiales;f__Bradyrhizobiaceae;g__Bradyrhizobium | 0.0064 | 0.0029 | 0.0078 | 0.0034 | 0.0035 | 0.0113 |
p__Proteobacteria;o__Enterobacteriales;f__Enterobacteriaceae;g__Escherichia | 0.0018 | 0.0104 | 0.0229 | 0.0055 | 0.0058 | 0.0071 |
p__Firmicutes;o__Clostridiales;f__Clostridiaceae;g__Clostridium | 0.0025 | 0.0247 | 0.0047 | 0.0009 | 0.0014 | 0.0007 |
p__Proteobacteria;o__Burkholderiales;f__Comamonadaceae;g__Variovorax | 0.0060 | 0.0047 | 0.0083 | 0.0029 | 0.0072 | 0.0017 |
p__Firmicutes;o__Clostridiales;f__;g__ | 0.0025 | 0.0247 | 0.0047 | 0.0020 | 0.0019 | 0.0009 |
p__Firmicutes;o__Turicibacterales;f__Turicibacteraceae;g__Turicibacter | 0.0009 | 0.0049 | 0.0016 | 0.0013 | 0.0002 | 0.0015 |
TS, piglets born to treated sows remained nontreated; TP, piglets born to non-treated sows were treated with ceftiofur; Control, non-treated piglets born to non-treated sows.
At D21, the number of differential ASVs was higher when comparing the control group with TS (92 ASVs) than with TP (68 ASVs). Among the ASVs affected by the antibiotic treatment, we found an ASV from Mycoplasma genus that was decreased in TS and absent in TP. This relative decrease in Mycoplasma in TP group was assigned, at least partially, to M. hyorhinis by PCR (Additional file 3: Table S2). Differentially abundant ASVs among the three groups (TS, TP and control) belonged to nasal-associated taxa, such as Streptococcus (more relatively abundant in TP), Moraxella and Glaesserella (both more relatively abundant in the control group). Interestingly, the same evidence was not observed at the genus level, suggesting that the antibiotic treatment in sows or piglets selected specific strains in each case. Within Glaesserella genus, G. parasuis is the only known member of the swine microbiota, for which a virulence-specific PCRs is available. We observed different colonization dynamics by virulent or non-virulent strains in the three groups, especially at early time points (Additional file 3: Table S2). The differences in ASVs were also supported by the isolation of colonies with a different fingerprinting profile in ERIC-PCR in the different groups. At genus level, not significant differences were detected, suggesting that the distribution of the different ASVs compensate each other yielding similar relative abundances in the groups. At D49, a global increase in the number and abundance of ASV was evident in all groups (as indicated in alpha diversity) and corresponded mostly to ASVs from Bacteroidales (9.40% in TS, 13.01% in TP and 12.56% in control) and Clostridiales (23.79% in TS, 21.54% in TP and 40.63% in control). At the end of nursery 593 ASVs were differentially abundant among the three groups (Additional file 4: Table S3). At genus level, some differentially abundant genera showed to be group-specific, such as Mycoplasma, whose abundance was higher in the control group (1.28% versus 0.15% in TP and 0.33% in TS). In agreement, higher prevalence of M. hyorhinis was detected in the control group by PCR (95% versus 74% in TP and 66% in TS). Sixteen genera were differentially abundant comparing the three groups (TS, TP and control), including some nasal-associated commensals: Actinobacillus (higher in TP), Bordetella (higher in TS), Lactobacillus (higher in TP), Staphylococcus (higher in TS) and a genus from the Pasteurellaceae family (higher in TS) (Additional file 5: Table S4). Ceftiofur treatment did not affect the prevalence of S. suis as detected by PCR, with a prevalence of over 90% at D7 and 100% at weaning and onwards in the three groups.
Inoculation of early colonizers of the upper respiratory tract modifies the nasal microbiota of the piglets.
Inoculation of non-treated piglets with selected colonizers (IP group) had an impact on their nasal microbiota. At D7, the mean alpha diversity measured by Shannon diversity index was higher in the inoculated IP piglets (Fig. 1) compared to the control group, although not significant at any timepoint (Kruskal-Wallis, P = 0.075, 0.14 and 0.25 for D7, D21 and D49, respectively). The microbiota composition estimated through beta diversity based on weighted Unifrac distance matrix was calculated for each timepoint (Fig. 4). Control and inoculated groups showed statistically different beta diversity at each timepoint (PERMANOVA, P = 0.02, P = 0.01 and P = 0.043 for D7, D21 and D49, respectively). The percentage of explanation attributed to the inoculation of the piglets (estimated through the Adonis function on the weighted Unifrac distance matrix) was 53% at D7, 44% at D21 and 36% at D49. Differential abundance analysis performed with ANCOM-BC methodology showed a total of 379 differential ASV, including four of the ASVs from the inoculated bacteria Moraxella, Rothia, Streptococcus and Glaesserella, which showed higher abundance in the IP group compared with the control group at D7 (Fig. 5), indicating that four out of five strains inoculated were able to colonize the nasal cavity of the piglets. The ASVs corresponding to the inoculated Glaesserella, Streptococcus and Moraxella were also present at weaning (D21), although significative differences were only detected with the Streptococcus ASV. These ASVs were not significantly different at the end of nursery where only Glaesserella and Streptococcus were detected at this latter timepoint. The inoculated Glaesserella strain was also recovered by culture at all the timepoints.
The colonization of the piglets by the selected colonizers seemed not to be favored by treating the sows with ceftiofur. When alpha diversity (measured through Shannon index) of inoculated piglets from treated or non-treated sows (IP vs TS-IP) was compared, no significant differences were detected at any of the time points analyzed (Kruskal Wallis, P > 0.05; Fig. 1). Differences in microbiota composition (beta diversity) were not detected with any of the metrics/distances analyzed (Jaccard, Bray Curtis and Unifrac; Fig. 6) at any of the sampling timepoints (IP vs TS-IP; Bray Curtis PERMANOVA, P = 0.396). The nasal microbiota composition of the two inoculated groups was not different at D7, although it was different for both groups when compared to the control (PERMANOVA, TS-IP vs control, P = 0.008 and IP vs control, P = 0.013, at D7). ANCOM-BC analysis showed a total of 379 ASVs at D7, including the ASV of the inoculated Rothia (0.0006% in IP and 0.0002% in TS-IP, while absent in control group). Similar findings were observed at D21 (PERMANOVA, IP vs TS-IP, P = 0.142; IP vs control P = 0.011; TS-IP vs control, P = 0.007) in beta diversity analysis, but at this timepoint none of the colonizers were differentially abundant. At the end of nursery (D49), inoculated groups were not significantly different (IP vs TS-IP; Bray Curtis PERMANOVA, P = 0.062). At this timepoint, only Glaesserella and Streptococcus were detected in the samples, but no significant differences were found in the differential abundance analysis (ANCOM-BC). Inoculation of the selected nasal colonizers reduced the divergence of the microbiota composition produced by the treatment of the sows. The differences detected at D7 and D21 in beta diversity were absent at the end of the study, indicating TS-IP piglets had a similar composition to the control (P = 0.163, with all beta diversity metrics used, while piglets born to treated sows but not inoculated (TS) presented different composition than control at D49. The IP group also showed a different microbiota than control at D49 (P = 0.008). Interestingly, the inoculated colonizers were found in higher relative abundance in IP than in TS-IP.
WGS metagenomics confirmed the differences observed in beta diversity between weaning and the end of nursery in all groups.
A total of 143 species and 99 genera were identified in the nasal microbiota samples from D21 and D49. In the control group, we identified a total of 10 genera and 14 species at weaning (D21), and the composition of the nasal microbiota was clearly modified through the incorporation of other members at D49 (Fig. 7A Top 50 species z-score control group D21 y D49). The divergence between D21 and D49 was also observed in piglets from the rest of the groups in a clustering analysis by average abundance (Z-score) at different taxonomic levels, including species (Additional file 6: Fig. S2, Bray Curtis top 25 species D21, D49 in Ts, TP and control). At D21, the top 10 most abundant genera in all groups were Glaesserella, Moraxella, Bergeyella, Mycoplasma, Streptococcus, Pasteurella, Neisseria, Mannheimia, Lactobacillus and Acinetobacter, which belong to taxa commonly associated to nasal microbiota. At the end of the nursery (D49), a plethora of species from taxa classically associated to gut microbiota and absent at weaning was detected. The same colonizers were detected at weaning as described for the 16S sequencing data and expanded their detection at the end of nursery detecting not only the Rothia at the end of nursery but Glaesserella strain as well.
Several species within Archaea domain were affected by the antibiotic treatment. While they were not found in any group at weaning, at the end of the nursery they were found more relatively abundant in the control and the IP group than the piglets born to treated sows (TS and TS-IP), although in less than one percent in all groups. Archaea were absent in the TP group at this latter timepoint as well.
A total of 314 pathways and 1580 enzymes were found considering all samples. Most of the pathways were present in all groups (Additional file 7: Fig. S3), except TS group which showed several specific pathways at D21, although represented in very low relative abundance. The differences found at taxa level between weaning and the end of nursery were also found in the functional profiling but were not so remarkable (Fig. 7B).
In the control group, 12 significantly different pathways were detected (through multivariable association analysis with linear model) when compared longitudinally (between D21 and D49 timepoints).These pathways were not among the top-50 most abundant pathways, and included several pathways in relative higher abundance at D21 (peptidoglycan maturation, heme b biosynthesis II, and two lipid IVA -precursor of lipid A- biosynthesis pathways) and others with higher abundance at D49 (anaerobic degradation of acetylene, degradation of stachyose, degradation of fucose and rhamnose, incomplete reductive TCA cycle, de novo biosynthesis of NAD from aspartate, L-arginine biosynthesis and two pathways of biosynthesis of pyrimidine).
Focusing on differences due to the antibiotic treatment at each time point, we detected few differential pathways, which were in higher abundance in the TS group (L-glutamate and L-glutamine biosynthesis at D21 but also at D49, superpathway of anaerobic sucrose degradation at D21, glycolysis-TCA-Glyox-bypass: superpathway of glycolysis at D49, and L-glutamine biosynthesis II at D21) when compared to TP and control groups.
In inoculated groups, we observed several pathways differentially increased at D21. Interestingly, these pathways increased in all groups at D49, indicating that these functions appeared earlier due to intervention, such as stachyose degradation, NAD de novo biosynthesis from aspartate, dTDP-beta-L-rhamnose biosynthesis (Additional file 8; Fig. S4). On the other hand, three pathways were found at lower abundance in both inoculated groups at D21 when compared to the control group, showing a delay in their natural appearance: superpathway of branched chain aminoacid, L-isoleucine biosynthesis I (from threonine) and L-isoleucine biosynthesis III. In addition, some differences were specifically observed in each inoculated group. The IP group showed two differentially abundant pathways: one in higher abundance at D49 (chitin derivatives degradation), and a second with lower abundance also at D49 (beta-(1,4)-mannan degradation). In TS-IP group, we observed five pathways that were more abundant at D21: incomplete reductive TCA cycle, NAD salvage pathway II (PNC IV cycle), 4-deoxy-L-threo-hex-4-enopyranuronate degradation, superpathway of arginine and polyamine biosynthesis, and purine nucleobases degradation II (anaerobic), where the latter three pathways appeared in lower relative abundance at D49 in the same group TS-IP (Additional file 8; Fig. S4).
The divergence observed in the taxonomic profiles among the groups was not so evident at the functional level, suggesting that different microbial communities may be responsible for similar functions.