3.1 Profiles of metagenome sequencing in the colonic digesta of growing pigs
A total of 6796902028 reads were generated via shotgun metagenome sequencing with 151042267 ± 2446487.48 (mean ± standard error of the mean) reads per sample (Additional file 2). After filtering the host genes and performing quality control, a total of 6498723288 clean reads were obtained with 144416073.1 ± 2417188.85 clean reads per sample. Then a total of 15115793 contigs (335907±11948 per sample) were obtained followed by de novo assembly.
3.2 Dynamic fluctuations of the colonic microbiota
These clean reads were further annotated with the NR database. After BLAST alignment against the NR database, these contigs were mainly annotated as bacteria (99.68 ± 0.0064%) or eukaryotes (0.23 ± 0.0070%) (Additional file 3). Using JTK_circle, a nonparameter analysis (Figure 1A), the relative abundances of bacteria (PAdj = 4.09 × 10-3) and eukaryotes (PAdj = 2.23 × 10-7) significantly fluctuated throughout 24h. Interestingly, bacteria were greater at night than in day in abundance (P = 0.003, Figure 1B). In contrast, that of eukaryota was boomed in day (P = 0.039, Figure 1B).
The microbial communities of bacteria and eukaryotes of different sampling timepoints exhibited obvious differences according to the PCoA model (Figure 1C). For bacteria, the dominant phyla were Firmicutes (80.64 ± 1.47%), Bacteroidetes (15.61 ± 1.40%), Proteobacteria (1.03 ± 0.057%), and Actinobacteria (0.96 ± 0.046%) (Additional file 3). Notebly, all of these parameters fluctuated significantly (PAdj < 0.05). Of note, the average cumulative relative abundance of these cycling phyla was over 98%. Among the top 10 most dominant genera (Figure 1D), the relative abundances of Lactobacillus (PAdj = 1.89×10-3), Prevotella (PAdj = 3.88 × 10-6), Bacteroides (PAdj = 3.42 × 10-5), Clostridium (PAdj = 0.021), Roseburia (PAdj = 3.35×10-3), Butyricicoccus (PAdj = 1.77 ×10-5), and Subdoligranulu (PAdj = 0.021) significantly fluctuated. At the species level, 19.29% of the annotated species exhibited significant fluctuations (PAdj = 0.05) in relative abundance (Figure 1E, F). Strikingly, the cumulative relative abundance of these cyclical species reached 58.86% (Figure 1F). Among the top dominant species, the relative abundances of Lactobacillus johnsonii (PAdj = 6.26 × 10-3), Lactobacillus reuteri (PAdj = 1.74 × 10-3), Prevotella copri (PAdj = 3.38 × 10-7), Butyricicoccus porcorum (PAdj = 5.56 × 10-6), Ruminococcaceae bacterium (PAdj = 2.85 × 10-3), Subdoligranulum sp. (PAdj = 9.11 × 10-3), and Gemmiger formicilis (PAdj = 6.76 × 10-3) exhibited significant oscillations (Figure 1E).
For eukaryotes, Fungi (27.42 ± 1.01%), Metazoa (41.02 ± 0.85%), and Viridiplantae (26.45 ± 0.88%), accounted for the largest proportion at the kingdom level (Additional file 3). Remarkably, all of them exhibited significantly oscillations. At the phylum level, except for Mucoromycota, the other 8 phyla Nematoda (PAdj = 2.89 × 10-4), Streptophyta (PAdj = 7.95 × 10-5), Ascomycota (PAdj = 3.49 × 10-4), Chytridiomycota (PAdj = 5.56 × 10-6), Chordata (PAdj = 7.29 × 10-3), Preaxostyla (PAdj = 0.021), Parabasalia (PAdj = 1.74 × 10-3), and Evosea (PAdj = 0.011) exhibited significant rhythmicity. In addition, 15.22% of the annotated species in eukaryotes underwent rhythmic oscillations (Figure 1E, F). The accumulated relative abundance of which contributed 39.74% to the total number of eukaryotes (Figure 1F).
Interestingly, rhythmic bacterial species peaking at different times within a day were mainly clustered into 7 classes (Figure 2A). Over 50% of Lachnospiraceae species were clustered into class 3 which had a trough during T15~T21, over 50% of Bacteroidaceae and Prevotellaceae were clustered into class 4, which had a peak during T06~T12; and over 50% of Erysipelotrichaceae were clustered into class 5, which had a peak during T12~T15. Over 50% of Clostridiaceae, Eubacteriaceae, Lactobacillaceae, and Ruminococcaceae were clustered into class 5, whose abundance peaked mainly during T15~T21 (Figure 2B).
3.3 Dynamic fluctuations in the amino acid profile of the colonic digesta
Data from untargeted metabolomics revealed that the metabolite profiles (including AAs, sugars, and fatty acids) in the digesta underwent rhythmic fluctuations over a day (Additional file 4). The metabolome explained the variation in microbial composition to a certain extent (Additional file 5). AA profile explained 13% of the variation in bacterial community composition, whereas organic acid profile explained 5% of the variation in eukaryotic community composition (Additional file 4, Additional file 5).
Therefore, we further determined AA profile in the colonic digesta by the targeted metabolomics. These results indicated that AA profile exhibited robust fluctuations within the day in the PCoA model. The T09, T12, T15, and T18 samples were distinguished from the T21, T24, and T27 samples (Figure 3A). All 35 identified amino acids and their derivatives except L-homoserine (PAdj = 0.553) exhibited significant fluctuations (Figure 3B). The most rhythmic AAs were citrulline (PAdj = 5.39×10-10) and L-methionine (PAdj = 9.62×10-10). Furthermore, most of these AAs were positively correlated with feed intakes (Figure 3C). These rhythmic AAs peaked at different times within a day and were mainly clustered into 5 classes (Figure 4A). The branched-chain amino acids L-valine (PAdj = 7.88 ×10-4), L-leucine (PAdj = 1.89×10-4) and L-isoleucine (PAdj = 1.85 ×10-4) were clustered into class 2, which peaked at T06~T12 (Figure 4A, Figure 5). The concentrations of L-valine (P = 0.011), L-leucine (P = 4.06×10-5) and L-isoleucine (PAdj = 3.02 ×10-5) were higher in dark phase than during light phase (Figure 4B, Figure 5). The aromatic amino acids L-tryptophan (PAdj = 2.26 ×10-3), L-tyrosine (PAdj = 4.50 ×10-6) and L-phenylalanine (PAdj = 8.61 ×10-7) were clustered into class 4 which were peaked at T12~T18 (Figure 4A, Figure 5). The most rhythmic AA citrulline (PAdj = 5.39×10-10) and L-methionine (PAdj = 9.62×10-10) were clustered into class 5 which were peaked at T09~T15 (Figure 4A, Figure 5). The concentrations of citrulline (P = 1.13×10-3) and L-methionine (P = 0.034) were higher in light phase than during dark phase (Figure 4B, Figure 5).
3.4 Amino acids dynamics affected the fluctuation of the gut microbes and the resistome
Based on a dynamic-correlation model eLSA, general relationships among the AA profile, the gut microbes, and the resistome were established. For the gut microbiota, most of these AAs positively correlated with Prevotellaceae species whereas negatively correlated with Lactobacillaceae and Lachnospiraceae species (Figure 6A). For instance (Figure 6B), citrulline (LS = 0.93, PAdj = 0.001, D = 2), 4-hydroxyproline (LS = 0.91, PAdj = 0.001, D = 3), and L-alanine (LS = 0.87, PAdj = 0.006, D = 2) were positively correlated with Prevotella copri whereas L-methionine (LS = -0.92, PAdj = 0.004, D = 0), and 4-hydroxyproline (LS = -0.86, PAdj = 0.012, D = 0) were negatively correlated with Lactobacillus reuteri, L-citrulline (LS = -0.91, PAdj = 0.005, D = -1), L-histidine (LS = -0.94, PAdj = 0.001, D = -1), and 4-hydroxyproline (LS = -0.89, PAdj = 0.007, D = 0) were negatively correlated with Roseburia sp. AM51-8.
The diurnal fluctuation in MGEs and ARGs were fully described in our work [15]. For the MGEs, AA profile negatively correlated with the Tn916 and istA group whereas positively correlated with plasmid (Figure 7A). For instance (Figure 7B), citrulline were positively correlated with rep14 (LS = 0.90, PAdj = 0.011, D = 1) and IncQ1 (LS = 0.83, PAdj = 0.021, D = 0). L-canisurine, an intermediate product of phenylalanine metabolism, were negatively correlated with Tn916-orf16 (LS = -0.97, PAdj < 0.001, D = 0), Tn916-orf18 (LS = -0.97, PAdj < 0.001, D = 0), and Tn916-orf20 (LS = -0.97, PAdj < 0.001, D = 0). For the ARGs, AAs had tightly correlations with ARGs belonging to major facilitator superfamily (MFS) antibiotic efflux pump family, tetracycline-resistant ribosomal protection protein, OXA beta-lactamase, and tetracycline-resistant ribosomal protection protein (Additional file 6). For example, citrulline were positively correlated with tet32 (LS = -0.89, PAdj = 0.007, D = -1) whereas negatively correlated with tetQ (LS = 0.90, PAdj = 0.004, D = -1) and tetB(P) (LS = 0.90, PAdj = 0.005, D = 1). L-canisurine were positively correlated with tetO (LS = 0.90, PAdj = 0.011, D = 1) and tetB(P) (LS = 0.93, PAdj = 0.002, D = 2) whereas negatively correlated with tetW (LS = -0.91, PAdj = 0.003, D = 0) and tetM (LS = -0.97, PAdj < 0.001, D = 0). Methionine was positively correlated with tetQ (LS = 0.93, PAdj < 0.001, D = -2) and tetB(P) (LS = 0.97, PAdj < 0.001, D = 1) whereas negatively correlated with tetM (LS = -0.95, PAdj < 0.001, D = -1) (Additional file 6).
Furthermore, using the partial least squares structural equation model (PLS-SEM, Figure 8), we found that AA profile directly affected the microbial α-diversity (R2 = 0.12), β-diversity (R2 = 0.60) and MGEs (R2 = 0.60). AA profile indirectly affected ARGs by affecting MGEs (R2 = 0.35) and microbial β-diversity (R2 = 0.35).