Understanding host-microbiota interactions in the commercial piglet around weaning

This work was part of a larger study that has been previously published^15,16,17 that is recommended for complementary information. In those publications we evaluated the impact of different management systems (including or not early socialization and environmental enrichment) on behavioural response¹⁵, performance¹⁶, and intestinal physiology and caecal microbiota¹⁷ of piglets. In the present work, our objective was to elucidate the changes induced by the commercial early weaning itself attending to changes on caecal microbiota, intestinal gene expression and metabolomic response. We also aimed to integrate all these data in a holistic approach to better understand the relationships between microbiota and animal response.

Weaning-induced changes in gut bacterial microbiome

An average of 78,562 ± 24,539 16S rRNA gene V3-V4 regions sequences per sample (ranging from 40,061 to 132,201) with an average length of 460 bp were obtained from the 28 caecal content samples, with no differences in the number of reads between pre- and post-weaning piglets (P = 0.424) and rarefaction curves reaching the plateau phase. The sequences were assigned to 976 Operational Taxonomic Units (OTU) based on a 97% sequence similarity. The indices of Chao1, observed species, Shannon and Simpson were calculated to estimate alpha diversity. As presented in Fig. 1., higher values were found after weaning (P = 0.015, P = 0.017, P = 0.013, P = 0.080; for Chao1, observed species, Shannon and Simpson indices, respectively), indicating an increase in the complexity of the gut ecosystem when animals are moved to the dry feed. Regarding beta diversity, a tendency was detected with the Whittaker index for a decrease after weaning (P = 0.062) indicating that ecosystems were more similar and a microbiota trend to converge between animals after weaning. Moreover, analysis of possible changes in the ecosystem structure related to weaning were performed using Anosim, Adonis and Envfit tests, all of them based on Bray–Curtis distance. Highly significant differences between pre- and post-weaning piglets were found (P = 0.0001, P = 0.001 and P = 0.0001, for Envfit, Anosim and Adonis tests, respectively).

Changes promoted in particular taxonomic groups

Figure 2 shows the relative abundances obtained for suckling and weaned piglets at phylum and genus levels. Firmicutes and Bacteroidetes constituted the two predominant phyla in the caecal microbiota of the piglets, followed by Proteobacteria (7.76%), Spirochaetes (3.49%) and Fusobacteria (2.85%). At genus level, Prevotella was found as the most predominant genus (13.61%), followed by unclassified Prevotella (7.70%) and Bacteroides (4.79%) from Bacteroidetes phyla.

Regarding the effect of weaning on microbial groups, although the most abundant phyla, Firmicutes, Bacteroidetes and Proteobacteria revealed no significant shifts in their relative abundances as a whole, significant effects were seen in predominant family and genera within. Only Fusobacteria phylum showed a remarkable decrease after weaning (P = 0.005). At the family level, four predominant families showed significant reductions as the piglets were weaned (Table 1), including Bacteroidaceae (P = 0.031), Enterobacteriaceae (P = 0.031), Fusobacteriaceae (P = 0.005), and Lactobacillaceae (P = 0.003). At the same time, Lachnospiraceae and Erysipelotrichaceae increased significantly after weaning (P = 0.031 and P = 0.022, respectively). At the genus level, two predominant genera showed significant increases after weaning (Fig. 3), including p-75-a5 (P = 0.001) and Roseburia (P = 0.012). The relative abundances of Bacteroides, Fusobacterium, Lactobacillus, and Megasphaera decreased in weaned piglets. Notably, Bacteroides and Fusobacterium, which were abundant in the gut of suckling piglets, declined from 7.01 and 4.72% to 2.58 and 0.99%, respectively for Bacteroides (P = 0.019) and Fusobacterium (P = 0.006), in a 5-day period. Among other non-predominant genera, several significant shifts were also detected, as for example, an increased abundance of Ruminococcus, Coprococcus, Dorea and Lachnospira in weaned piglets (P = 0.020; P = 0.018; P = 0.038; and P = 0.002, respectively).

Predicted functions of the caecal microbiota

The inference of the functional profile of the caecal microbial community was predicted by using PICRUSt¹⁸ v1.1.3. A clear differentiation was observed between lactating and weaned piglets’ microbiota related to several KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways¹⁹. During lactation, pathways related to metabolism, like glycan biosynthesis or lipid metabolism, cofactor and vitamin and nucleotide metabolism, were highly represented. After weaning, functions related to cellular processes, as cell motility and sporulation, or related to environment information processing, like membrane transport and signal transduction were higher (Fig. 4).

Going to a deeper level (KEGG level 3, presented in Supplementary Fig. S1), lipid metabolism pathways, such as lipid biosynthesis proteins (P = 0.005), energy metabolism pathways, such as carbon fixation pathways (P = 0.019) and carbohydrate metabolism pathways, such as the citrate cycle (TCA cycle, P = 0.028), among others, were more represented in suckling piglets. The pathways involved in metabolism of purine (P = 0.041) and alanine, aspartate and glutamate metabolism (P = 0.049), related to nucleotide and amino acid metabolism pathways, respectively, were also higher during lactation compared to suckling. Nicotinate and nicotinamide metabolism (P = 0.008), related to the metabolism of cofactors and vitamins, increased also in weaned piglets compared to lactation.

Weaned piglets showed, however, a higher proportion of pathways involved in bacterial chemotaxis (P = 0.004), bacterial motility proteins (P = 0.007) and flagellar assembly (P = 0.014), all related to cell motility. Pathways related with membrane transport, such as transporters (P = 0.018) and ABC transporters (P = 0.015), and signal transduction, such as the two-component system (P = 0.015) were higher after weaning. Finally, lipopolysaccharide biosynthesis proteins (P = 0.040), related to the glycan biosynthesis pathway, and sporulation pathway (P = 0.014), showed higher values in weaned piglets.

Changes induced in the jejunal gene expression

Jejunum samples from the fourteen piglets were collected shortly before and after weaning to analyse the expression of several genes related to intestinal functionality by using the OpenArray technology. The 51 genes analysed were grouped into six categories for easier understanding according to whether they were related to: barrier function (BF), immune response (IR), nutrient transport (NT), enzyme/hormone encoders (EH), stress indicators (ST) or housekeeping (HK).

Several genes from all functional groups showed significant changes after weaning as shown in Table 2. Moreover, some additional genes related to barrier function, tended to be downregulated (CLDN4 and MUC2, P = 0.064 and P = 0.074 respectively) or increased (CLDN15, P = 0.079) in weaned piglets.

Weaning-induced changes in serum metabolome

¹H-NMR spectra

Twenty-eight serum samples were prepared and analysed from 14 suckling piglets and 14 weaned piglets. Among the different endogenous metabolites assigned there were LDL, VLDL, lipids, unsaturated lipids, leucine, valine, isoleucine, lactate, alanine, adipate, acetate, N-acetyl glycoproteins, O-acetyl glycoproteins, glutamine, glutamate, pyruvate, creatine, choline, trimethylamine-N-oxide (TMAO), glucose, creatinine, tyrosine and phenylalanine.

To identify potential differences between serum metabolites profiles of pre- and post-weaning piglets, an untargeted metabolomics approach using ¹H-NMR was also applied. In order to reduce the number of variables, a filtering of ¹H-NMR bucket table was done by significant differences on Student’s t-test (P-value ≤ 0.2) between the integrated buck regions of suckling and weaned piglets. To evaluate the global metabolic profile of serum samples collected from piglets in both periods, a blinded to age groups study by principal component analysis (PCA) of ¹H-NMR datasets was performed from the filtered ¹H-NMR bucket table. Figure 5a shows a biplot representation of PCA [R²x_(cum)=0.82, Q²_(cum)=0.30] from the reduced data, in which a clear pattern of separation between suckling and weaned piglets along PC1 could be observed, indicating that both piglets’ groups were metabolically differentiated. To identify the key metabolites that influence in this grouping, a study taking account the age groups was made by an orthogonal partial least squares discriminant analysis (OPLS-DA) approach (Fig. 5b). The supervised OPLS-DA model [R²x_(cum) = 0.38, R²y_(cum) = 0.84, Q²_(cum) = 0.60] developed a perfect separation into the two clusters with high fitness R² and accepted predictive ability Q² parameters (R²y_(cum) and Q²_(cum) > 0.5). Moreover, the cross-model validation (Supplementary Fig. S2a) and the permutation test (100 times) (Supplementary Fig. S2c), both indicate that the developed OPLS-DA approach was positive and valid, confirming the distinction among both piglets’ groups. Also, a value of 1.0 for the area under the curve (AUC) corresponding to receiver operating characteristic (ROC) plot (Supplementary Fig. S3) indicated a strong discrimination power for the OPLS-DA classifier model. To find the most relevant ¹H-NMR regions that contribute to the differentiation between suckling piglets from weaned piglets, an S-plot was performed (Supplementary Fig. S4). From this plot, the key ¹H-NMR buckets that affect the discrimination were identified. These regions were also screened according to their corresponding variable importance in the projection (VIP) values of the OPLS-DA model. Fifteen from the total 250 spectral regions were found like the more contributing (Table 3), of which 12 regions had integral values that differed significantly between both groups (Student’s t-test P-value ≤ 0.05). The discriminant metabolites that showed higher levels in suckling piglets were choline, lipids (including triglycerides and fatty acids), LDL, alanine, isoleucine and probably also TMAO, whereas in weaned piglets those with higher levels were 3-hydroxybutyrate, ethanol, valine, and adipate (Table 3).

Integration of the omics technologies

Gene expression, caecal microbiota and metabolomics datasets were integrated by using the open-source software R²¹ v3.6.1. and the LinkHD²² package. The objective of this approach was to analyse these heterogeneous datasets to verify from a holistic point view that weaning was determinant defining different clusters of samples. Therefore, confirming our hypothesis, with the additional value to explore the connections (i.e.: correlations) between the gene expression, metabolomics and cecal microbial communities. Furthermore, the use of LinkHD allow us to highly the most informative variables within each dataset, as well as to identify which dataset was most relevant for the sample stratification. Although no model was applied in the statistical design, the samples were stratified into two clusters in a blind analysis, that aligned quite well with the groups of suckling and weaned piglets. The different ordination was mainly explained by the changes observed in piglet microbiota but also by the differential distribution of metabolites. Gene expression did not appear to contribute significantly for the cluster ordination probably due to the low number of input variables (52 genes) compared with the microbial and metabolites datasets. In total, 93 OTUs and 12 metabolites were found as discriminating between groups although the metabolites were not able to be identified. The discriminant relevant families included among others, Fusobacteriaceae (P = 0.0010), Bacteroidaceae (P = 0.0011), Enterobacteriaceae (P = 0.0012), Lactobacillaceae (P = 0.0013), Erysipelotrichaceae (P = 0.0138) and Prevotellaceae (P = 0.0424). Meanwhile, at genera level the significant discriminant genera were among others, Fusobacterium (P = 0.0010), Bacteroides (P = 0.0011), Lactobacillus (P = 0.0013), Megasphaera (P = 0.0014) and Ruminococcus (P = 0.0157).

In addition to the stratification in clusters and the identification of the most relevant variables, LinkHD also related the existing multivariate correlation (RV) between the different datasets (Fig. 6). Thus, it was observed that the strongest association was that between the metabolites and 16S (RV values: Genes & 16S = 0.37; Genes & Metabolites = 0.22; 16S & Metabolites = 0.53).

Integration of gut microbiome and serum metabolome data

Considering the relevance of microbiota and metabolomic data in the LinkHD analysis, an additional integration approach was performed by correlating the caecal microbiota with the blood serum metabolome. The objective was to elucidate possible relationships between the intestinal microbiota and the animal metabolism, exploiting the wide range of response found around weaning. A multivariate analysis of the data was performed looking for significant correlations among the traits (Supplementary Table S2). The variables integrated were ¹H-NMR bucket area regions of 0.04 ppm width (58) and bacteria family counts (22), obtained from nursing piglets (n = 14) and weaned piglets (n = 14). The 58 bucket regions, from the total of 250 regions of the spectra, were selected based on previous assigned soluble metabolites by Clausen et al.²³ and He et al.^24,25 The metabolites that showed significant Pearson correlation coefficients |r|≥ 0.37 (P-value ≤ 0.05) with the relative abundance of bacterial at family taxonomic level can be found in Supplementary Table S3.

Furthermore, to better understand the putative relationships between gut microbiota and serum metabolites, a global analysis was performed by a hierarchical clustering of Pearson correlations coefficients (Fig. 7). This analysis evidenced four clusters for metabolite buckets (illustrated as cluster I-IV in Fig. 7) and 3 major clusters for bacteria families (illustrated as A, B and C). Interestingly, the clusters evidenced for metabolite buckets corresponded to certain chemical or metabolic categories. In this way, Cluster I included mainly glucose; Cluster II appears to be related to protein metabolism including aromatic (Tyr, Phe) and other amino acids (Ala, Gln, Glu) and glucose metabolism (glucose, pyruvate and lactate); Cluster III mostly included branched chain amino acids (Ile, Val) and also compounds related to energy metabolism (creatin/creatinine, 3-hydroxybutyrate and acetate) and Cluster IV included buckets mostly related to lipid metabolism (unsaturated lipids, lipids [triglycerides and fatty acids], LDL, VLDL, choline). Regarding microbial clusters, Cluster A conformed by Mogibacteriaceae, Coriobacteriaceae, Lachnospiraceae, Clostridiaceae, Ruminococcaceae, Spirochaetaceae and Streptococcaceae families, showed to be in general terms positively correlated with Cluster I and cluster III and negatively correlated with Cluster IV. Cluster B, including Alcaligenaceae, Campylobacteraceae, Lactobacillaceae, Paraprevotellaceae, Porphyromonadaceae, Erysipelotrichaceae, Desulfovibrionaceae, S24-7, Prevotellaceae and Veillonellaceae, differed from Cluster A in that Cluster B was negatively correlated with Cluster I and positively correlated with Cluster II. Finally, Cluster C, including Bacteroidaceae, Dethiosulfovibrionaceae, Odoribacteraceae, Enterobacteriaceae and Victivallaceae families, differed from Clusters A and B in that it was strongly negative correlated with cluster I, II and III, and weakly positive correlated with cluster IV.