Nutritional analysis of ensiled and non-ensiled GM
The nutritional traits of GM were not altered with the ensiling processes but some numerical changes were rather observed with the treatments (Table 1). Crude proteins increased with the treatment D1+i and D2+i. No matter the initial DM (D1 or D2), the protein content increased when adding the inoculant compared to when adding the sucrose. DM immediately after ensiling did not change compared to the respective initial DM. The treatment D1+i had the highest protein content and was therefore selected for further metabolomic and sequencing analysis. pH value decreased and was finally maintained within the normal range of acidity for ensiling. Except for D1+S that increased the temperature, the other treatments were maintained. To further understand the effect of the inoculant on ensiling, sample 1, sample 2 and sample 3 standing respectively for ensiled treatment (D1+i), non-ensiled treatment and inoculant were investigated for the metabolomic analysis.
Table 1: Nutritional traits of ensiled and non-ensiled GM
Trt
|
CP (%)
|
NDF(%)
|
ADF(%)
|
Ash(%)
|
DM*(%)
|
DM#(%)
|
pH
|
Temp. (°C)
|
NE GM
|
13.89±0.6bc
|
39.81±3.0bc
|
24.46±1.5bc
|
13.98±0.49a
|
NA
|
94.58±0.62a
|
7.75±0.35a
|
NA
|
D1
|
14.06±0.19b
|
43.62±3.4a
|
27.11±0.72a
|
13.03±0.02b
|
50.87±2.01b
|
90.87±0.46bc
|
4.14±0.02c
|
20.0±1.35b
|
D1+i
|
15.43±0.22a
|
40.65± 3.2b
|
25.73± 0.92b
|
13.38±0.13b
|
51.07±1.73b
|
91.38±0.83b
|
4.28±0.03c
|
19.9±1.35b
|
D1+S
|
14.6±0.73b
|
38.66±3.1c
|
24.52±1.51bc
|
13.20±0.38b
|
50.58±1.93b
|
92.13±1.01ab
|
4.28±0.03c
|
21.6±1.21a
|
D1+i+S
|
14.12±0.18b
|
39.60±4.9b
|
25.95±2.91b
|
13.25±0.07b
|
50.89±1.88b
|
91.51±076b
|
4.26±0.01c
|
19.9±1.14b
|
D2
|
13.46±0.56c
|
36.58±2.7c
|
23.89±2.76c
|
13.60±0.12ab
|
59.87±0.45a
|
88.46±1.01c
|
5.38±0.41b
|
20.6±1.41b
|
D2+i
|
15.13±0.17a
|
/
|
/
|
13.33±0.49b
|
61.17±0.53a
|
89.88±0.26c
|
5.35±0.14b
|
20.4±1.01b
|
D2+S
|
13.92±0.43bc
|
41.76±4.4ab
|
26.65±1.12ab
|
12.87±0.16bc
|
61.00±0.49a
|
89.60±0.44c
|
5.26±0.35b
|
19.7±1.02b
|
D2+i+S
|
14.02±1.0b
|
33.01±0.7d
|
21.17±1.56d
|
12.67±0.67c
|
60.77±0.51a
|
88.65±1.27c
|
5.41±0.21b
|
19.9±1.03b
|
Trt: treatment, CP: crude protein, NDF: neutral detergent fiber, ADF: acid detergent fiber, Temp.: temperature, DM: dry matter, *: DM immediately after 60 days ensiling, #: DM after 48 h oven drying at 65° C and 4 h at 105° C, NE GM: non-ensiled GM, D1: DM 50 %, D2: DM 60 %, D1+I: DM 50 %+inoculant, D2+I: DM 60 %+inoculant D1+I+S: DM 50 %+inoculant+sucrose, D2+I+S: DM 60 %+inoculant+sucrose, D2+S: DM 60 %+sucrose, D1+S: DM 50 %+sucrose. NA: Not applicable, Different superscripts within columns denote a significant difference, P˂0.05.
Metabolomic profile of GM, ensiled GM and inoculant bacteria
The flavonoids, lipids, organic acid, nucleotides and derivatives content decreased in the ensiled GM compared to the non-ensiled GM. Lignans and coumarins, terpenoids and steroids were balanced while phenolic acids, amino acids and derivatives, alkaloids and quinones increased (Fig. 1a). The inoculant bacteria metabolites compared to GM, recorded a lower
X-axis indicates the sample name and the Y-axis are the class of metabolites. Sample 1, 2 and 3 are respectively the ensiled GM, the non-ensiled GM and the inoculant bacteria with 3 replications. Group indicates sample groups. Z-Score indicates the relative quantification of each metabolite with red representing higher content and green representing lower content.
content of flavonoids, lipids, phenolic acids, lignans and coumarins, terpenoids, steroids and quinones but a higher amino acid and derivatives compared to the non-ensiled GM (Fig. 1b). A detailed hierarchical cluster analysis of the metabolites at a lower level, before and after GM ensiling is shown in Figure S1.
K-Means and Venn diagram analysis of differential metabolites in the treatment groups
K-Means analysis (Fig. 2) examined the trend of relative quantification changes of metabolite sub-classes with different treatments or sample groups. A total of 10 sub-classes based on metabolic trends across the treatments were represented. Sub-class 1 recorded 77 metabolites with the same tendancy along with sample 1 and 2 but the highest relative quantification or abundance in sample 3. Sub-classes 3, 4, 5, 8 and 9 recorded respectively 78, 69, 100, 69 and 64 metabolites with the lowest abundance in sample 3. Sub-class 2 with 282 metabolites and sub-class 5 with 100 metabolites have the highest abundance in sample 2. After fermentation (sample 1), metabolites from sub-classes 3, 4, 6 and 8 have been found with the highest abundance.
The X-axis represents the sample names and the Y-axis represents the normalized relative quantification. ”Sub Class” represents a group of metabolites with the same trend and the number represent the number of metabolites in this cluster. K-Means is performed based on the Z-score normalized relative quantification value.
Venn diagram revealed that a major portion of the metabolites were overlapped among the different treatments, whereas few metabolites were stage-specific (Fig. 3). For instance, only 33, 60 and 23 metabolites were specific to samples 3, 2 and 1 respectively. The other inter-connected metabolites were detected among all the different samples. GM, ensiled GM and inoculant display non negligeable metabolic differences that remained to be characterized at a lower level.
Principal component analysis of the treatment groups
The principal component analysis (PCA) presented three distinct clusters (inoculant bacteria, non-ensiling plant and QC samples) observed along PC1, which represented 52.94% of the variation among the samples (Fig. 4), and suggested the level of metabolome change during ensiling. The ensiling plant and the inoculant bacteria were separated along PC2 with 38.09% of the variation. The separation between the ensiling and non-ensiling plant was not observable in this 2-dimensional PCA, neither along PC1 nor along PC2.
Dynamic distribution of the top 10 metabolites content
The top 10 up-regulated and down-regulated expressed metabolites were distributed in 3 classes or categories of fold-change log2FC values (Fig. 5). 3,4'-dihydroxybenzoic acid ethyl ester, 2-hydroxyethylphosphonic acid, 3,4'-dihydroxy-3',5'-dimethoxypropiophenone, vnilloylmalic acid constitute the first category with an up-regulation rate ranging from 13 to 14 log2FC. Sedoheptulose, 2-hydroxy-3,5-dinitrobenzoic acid and L-arginine were expressed at a rate of 16 to 17 log2FC while putrescine, methyl linolenate and calactin were expressed at the highest rate of 18 to 19 log2FC. Meanwhile, among the top 10 down-regulated expressed metabolites, 2’-o-methyladenosine was the most down-regulated metabolite with -19.5 log2FC. Xanthosine, 2-hydroxy-2-methyl propyl glucosinolate and isopentenyl adenine-7-N-glucoside have a similar down-regulated log2FC value around -17. The other down-regulated metabolites were around -16 (Fig. 5). The chemical formula, class and CAS numbers of the top 10 differential metabolites are presented in Table S1.
X-axis refers to log2FC values of top differential metabolites, the Y-axis refers to metabolites. Red bars represent up-regulated differential metabolites and green bars represent down-regulated differential metabolites.
Spearman correlation analysis
The inter-relationship between metabolites has been studied through a Spearman rank’s correlation (Fig. 6) and a Pearson’s correlation (Fig. S3) that revealed a synergistic or antagonistic relationship among the metabolites. The negative relationship between the down-regulated and the up-regulated metabolites on the one hand, and a positive relationship either between the up-regulated metabolites or between the down -regulated metabolites were observed. When the metabolism of 3,4'-dihydroxy-3',5'-dimethoxypropiophenone increased, the one of xanthosine decreased while L-arginine rather increased.
KEGG enrichment analysis of differential metabolites
Biological pathways and functions are made possible with the contribution of a set of metabolite populations that vary with the pathways (Fig. 7). Therefore, the highest number of metabolites (70), were involved in the "biosynthesis of secondary metabolites" pathway. ABC transporters, biosynthesis of cofactors and biosynthesis of amino acids have the same number of metabolite contributions (30). A small number of metabolites (5) contributed to the thiamine metabolism, sphingolipid metabolism, isoquinoline alkaloid biosynthesis and carbon fixation in photosynthetic organisms. The distribution of metabolites inducing different KEGG pathways related to metabolism, genetic information processing and environmental information processing is detailed in Figure S2.
Microbiota composition and structure in ensiled and non-ensiled GM plant
The bacterial profiles differed to a greater extent at lower (e.g., specie) than at higher (e.g., class) taxonomic ranks with a decrease of the relative abundance at a lower taxonomic rank (Fig. 8). The relative abundance of the majority of identified bacteria increased with ensiling. Indeed, at the class level, Bacteriodia, Clostridia, Negativicutes, Coriobacteriia, Anaerolineae and others increased with ensiling while Cyanobacteria and Alphaproteobacteria decreased. At the family level, Bacteroidaceae, Ruminococcaceae, Lachnospiraceae, Rikenellaceae, Prevotellaceae, Oxillospiraceae, Acidaminococcaceae and others increased in ensiled samples but Rhizobiaceae and Sphingomonadaceae decreased. At the genus level, Bacteroides, Faecalibacterium, Rikenellaceae_RC9_gut_group, Phascolarctobacterium, Megamonas, Prevotella increased with ensiling while Aureimonas and Pseudomonas decreased. At the specie rank, Bacteriodes salanitronis, B. plebeius, B.
barnesiae, B. vulgatus, B. caecicola, Prevotella copri, Megamonas hypermegale, Olsenella sp. and others increased with ensiling processes. Sphingomonas paucimobilis was found only in the non-ensiled plant while Lactobacillus buchneri was specific to the ensiled plant. The ensiled samples were free of fungi and mold populations meanwhile, a diversity was observed in non-ensiled samples with Lasiodiplodia margaritacea the most represented specie belonging to the Batryosphaeriaceae family and Dothideomycetes class (Fig. 9).
Only the top 10 taxons of abundance are displayed on the bar chart.The remaining proportions to get 100% on each bar were unclassified taxons.
Only the top 10 taxons of abundance are displayed on the bar chart.The remaining proportions to get 100% on each bar were unclassified taxons.
Phylogenetic analysis of the ensiled and non-ensiled microbiota
To further understand the relationships among the bacteria, a phylogenetic tree was constructed using representative sequences of bacterial genera (Fig. 10a) and fungi and mold genera (Fig. 10b). The tree enabled to observe the phylogenetic distance between the studied microbiota. Among the Firmicutes, four sub- groups of bacteria were genetically individualized with UCG-004 their common ancestor. Individual like Shuttleworthia was similar with agathobacter than with Erysipelotrichaceae_UGG-003.
Within the Proteobacteria group, Pseudomonas is the parent of all the members but is also very close to the member of another group, Desulfobacteriota (Desulfovibrio). The members of Bacteroidota can be separated into three sub-groups of similar genotypes. Methanosaeta had the highest phylogenetic distance among all the studied bacteria. The fungi and mold tree was predominated by the Ascomycota members, with Tremateia, Kodamaea, Metschnikowia, and Candida genetically very far from the other Ascomycota but closer to the Basidiomyocota members. Mucoromycota was the least represented group.
Functional prediction
The relative abundance of KEGG pathways was predicted using PICRUSt2, based on representative sequences and categorized according to groups of functional roles (Fig. 11). The more expressed functional group was related to the "brite hierarchies" while the least one was related to "human diseases and organismal systems". The relative abundance of the pathways related to "signal transduction and membrane transport" increased with the non-ensiling while "translation" and "replication and repair" increased with ensiling. Within the human diseases functional group, the ensiling process reduces the relative abundance of "infectious viral diseases", "neurodegenerative deseases" and "cancer overview". Moreover, within the metabolism group, the relative abundances of "xenobiotics biodegradation and metabolism", "metabolism of cofactors and vitamins", "energy metabolism" and "metabolism of other amino-acids" increased in the non-ensiled plant. "Nucleotide metabolism", "glycan biosynthesis and metabolism" and "carbohydrate metabolism" were increased with ensiling.
Cause and effect relationships: Network correlation between the microbiota and the differential plant metabolites
A network correlation analysis enabled to characterize the interaction between the most changing metabolites and the microbiota profiles (Fig. 12). Sphingomonas paucimobilis is positively correlated to isopentenyl adenine-7-N-glucoside and 5-hydroxymethylfurfural. Bacteriodes vulgatus was positively correlated to L-arginine and 3,4'-dihydroxy-3',5'-dimethoxypropiophenone and negatively correlated to medioresinol-4'-O-(6''''-acetyl)glucoside, 8-(benzyloxy)-1-naphthol, 2-(2-phenylethyl)chromones and xanthosine. Bacteroides plebeuis and Prevotella copri are positively correlated to L-arginine, 3,4'-dihydroxy-3',5'-dimethoxypropiophenone and vnilloylmalic acid. The dynamic of most of the metabolites was explained by the contribution of Lactobacillus buchneri, Bacteriodes salanitronis, Bacteriodes barnesiae, Bacteriodes caecicola and Olsenella sp. that were correlated to most of the studied metabolites and could serve in their synthetic forms as inoculant to manipulate the metabolic profile of GM plant. Fungi and mold contributed very little to the metabolites dynamic, except for Filobasidium oeirense and Rhizopus arrhizus which both contributed to up and down-regulation.
In the network diagram, bacteria are represented by turquoise circles, fungi are represented by bright coral circles, and metabolites are represented by cornflower blue circles; positive correlations are represented by red solid lines, and negative correlations are represented by gray solid lines; the size of the circle indicates connectivity, and the thickness of the line indicates correlation. The correlation filtering threshold: |cor|>0.5, pvalue<0.05