Fermentation and nutritional quality profiling of alfalfa silage
The ensiling characteristics, including pH and organic acid content, are presented in Figure 1. Both the LN (L. plantarum NF92) and LB (L. plantarum B90) treatments significantly reduced the pH (LN: 4.18; LB: 4.22) of the alfalfa silage compared with that of the CK group (4.82, Fig. 1A). In comparison, the lactic acid (LA) content of the LN- and LB-treated silages was greater (p < 0.01) than that of the CK (Fig. 1B), while the acetic acid (AA) content increased only in the LB-treated silages. The ratio of LA to AA had the same tendency as that of the LA content (Fig. 1D). In the present study, considerable dry matter (DM) loss (6.24%) occurred in the CK group (Fig. 1F). However, the LN and LB treatments obviously decreased (p < 0.0001) the loss of DM in the silage, especially in the LB treatment. The water soluble carbohydrate (WSC) content of the LN group was much greater than that of the CK group, while that of the LB group did not significantly change (Fig. 1G). In addition, the acid detergent fiber (ADF) content of the LN-treated silages was significantly lower than that of the fresh forage (FR) and CK silages (Fig. 1H). The neutral detergent fiber (NDF) levels of the LN- and LB-treated silages were significantly lower than those of the FR- and CK-treated silages (Fig. 1I). In conclusion, the LN and LB treatments improved the silage fermentation quality by decreasing the pH and increasing the LA content and improved the nutritional quality by decreasing the DM loss and indigestible fiber content.
The nitrogen composition of alfalfa silage reveals variability in proteolysis between treatments
Although the CP content of all treated silages remained comparable to that of FR (Fig. 2A), distinct variations were observed in the protein fractions compared to the composition of FR. After ensiling, the level of PA2 (Fig. 2B) in the silages was significantly greater than that in FR, while that in PB1 (Fig. 2C), PB2 (Fig. 2D), and PC (Fig. 2E) was notably lower. Compared with those of the control silage, both the LB- and LN-treated silages exhibited a significant increase of more than 20% in the PB1 level (Fig. 2C). No significant difference was observed in the contents of PB2 and PC between the CK silage and silages treated with LB and LN (Fig. 2D, E). After ensiling, the NH3-N level in the CK silage significantly increased, whereas it decreased in the LN- and LB-treated silages (p < 0.0001) (Fig. 2F). The amino acid-nitrogen (AA-N) and peptide-N levels in the silages were also significantly (p < 0.001) greater than those in the FR treatment (Fig. 2G, H). Compared to the CK silages, the LN and LB treatments significantly (p < 0.0001) decreased the AA-N level but increased the peptide-N level. The results suggest that the LN and LB treatments could improve the utilization of silage N by increasing the true protein (PB1) and peptide-N contents. Additionally, these treatments contribute to mitigating nutrient loss by reducing NH3-N contents. Nevertheless, the protease PrtR had no obvious influence on the protein fraction.
Changes in the protein abundance of alfalfa silage reflected by the proteome
To obtain a representative overview of the protein alterations, a fingerprint analysis of the proteins was performed by SDS-PAGE (Fig. 3A). Compared to those of the FR samples, the profiles of the ensiled samples were modified to different intensities. After ensiling, the intensities of proteins with molecular weights (MWs) in the ranges of 45-55 kDa and 70-170 kDa decreased, whereas those of 55-70 kDa and 15-25 kDa increased (Fig. S1). The intensity of proteins at 25-35 kDa did not significantly vary. Minor differences appeared between the CK silages and the LN- and LB-treated silages. To further detect the differentially expressed proteins (DEPs) between silages, in-depth TMT-based proteome characterization was performed. Our workflow identified and quantified 31049 unique peptide sequences that were assigned to 5109 proteins (Table S1). Principal component analysis (PCA) revealed distinct proteomes between the samples, with the first component accounting for 46.95% of the variance, distinctly separating the FR from the silages (Fig. 3B). The second component further delineated LN- and LB-treated silages from CK silage, explaining an additional 23.06% of the variance. However, the LN- and LB-treated silages clustered closely, exhibiting greater similarity to each other. The most abundant protein in the FR group had an MW ranging from 45 to 55 kDa, as detected by proteomic analysis (Fig. 3C). Notably, the predominant protein within this range was identified as ribulose bisphosphate carboxylase large chain (RuBisCO large chain, S4T017, and A0A072TKG8), constituting 62.02% of the total abundance. The increase in proteins with MWs in the range of 55-70 kDa in the CK group is shown in Figure 3D, among which proteins with protease activity, such as A0A072VYH6 (fungal proteinase A), A0A072TLC8 (carboxypeptidase), G7JXM9 (eukaryotic aspartyl protease family protein), and A0A396HXU8 (putative tripeptidyl-peptidase II), were upregulated. To capture disparities in protein abundance across different groups, the protein intensity of the proteome was extracted. The FR group had the highest total protein abundance (Fig. 3E, F). Among the top 10 most abundant proteins, RuBisCO large chain (S4T017) and small subunit (A0A072U6Q1) collectively constituted 24.91% of the total protein mass in the FR group. Notably, the RuBisCO large chain was strongly downregulated in the ensilaged group, particularly in the CK group (Fig. 3E). Pairwise comparisons were performed to filter potential DEPs, encompassing both upregulated and downregulated candidates. Compared with those in the FR group, 352 DEPs were downregulated in the CK group, accounting for 38.98% of the total FR mass, while 429 upregulated DEPs accounted for 9.47% of the total FR mass (Fig. 3G, J). Compared to those in the CK group, 167 DEPs (accounting for 3.96% of the total mass) were downregulated in the LN group (Fig. 3H, J), while 180 DEPs (accounting for 4.32% of the total mass) were downregulated in the LB group (Fig. 3I, J). These results underscore that ensiling affected the proteome of the samples, with additional modulations observed in the LB and LN groups, further influencing the silage proteome.
Distinct proteome features of alfalfa silage between treatments
To determine the differences in the proteome between CK silage and silages treated with LN and LB, a Venn diagram was generated. Within the subset of downregulated proteins, a single protein was shared across all groups. A total of 146 DEPs were common to both the LN vs. CK group and the LB vs. CK group. There were 351, 20, and 33 unique DEPs exclusive to the CK vs. FR group, LN vs. CK group, and LB vs. CK group, respectively (Fig. 4A). Among the upregulated DEPs, 427, 5, and 3 unique DEPs were distributed in the CK vs. FR group, LN vs. CK group, and LB vs. CK group, respectively (Fig. 4B). Given that the variation in the number of DEPs may not reflect their variation in abundance, the abundance of DEPs was determined. In the CK vs. FR group, DEPs with MWs of 45-55 kDa had the highest abundance (48.84%, 103 proteins), among which downregulated proteins accounted for 45.89% (Fig. 4C, F). In the LN vs. CK group (27.87%) and LB vs. CK group (22.81%), DEPs with MWs of 35-45 kDa had the highest abundance (Fig. 4D, E, G). Our data indicated that the ensiling process mainly downregulated the abundance of DEPs with MWs of 45-55 kDa. Additionally, the LN and LB treatments further downregulated the abundance of the 35-45 kDa DEPs in conjunction with the ensiling process.
Proteomic alterations are associated with protein cellular localization
To investigate the DEP composition of subcellular compartments, subcellular maps were generated according to Gene Ontology Cellular Component annotations. Many proteins lack subcellular localization annotation, and some proteins co-localize to more than one organelle (Fig. 5A). Membrane proteins comprised 40.35% of the total protein quantity, followed by cytoplasmic proteins (29.65%) and ribosomal proteins (13.63%) (Fig. 5B). To further dissect the subcellular localization of DEPs regulated by different treatments, the localization of DEPs in the CK vs. FR group, LB vs. CK group and LN vs. CK group was determined. Compared to those in FR silage, the DEPs whose levels were upregulated in CK silage were mainly located in the membrane (45.0%) and extracellular regions (40.0%), while those whose levels were downregulated were located in the ribosome (36.51%), membrane (30.95%) and cytoplasm (21.43%) (Fig. 5C, D). Given that the differences in the number of DEPs across subcellular locations may not reflect their abundance variation, the protein abundance was estimated (Fig. 5G). The abundance of membrane proteins regulated by ensiling accounted for the greatest proportion (78.37%), among which the downregulated proteins accounted for 44.40%. The abundances of ribosome and cytoplasmic proteins were 13.86% and 7.52%, respectively, while those of mitochondria and chloroplasts were the lowest. Compared to those in the CK group, the DEPs in the cytoplasm and membrane in the LN and LB groups were further downregulated (Fig. 5E, F). The abundance of cytoplasmic protein accounted for the greatest proportion, with 72.62% in both groups (Fig. 5H, I). This result indicated that the ensiling process primarily induced changes in DEPs located in the membrane, ribosome, and cytoplasm, while LN and LB further downregulated cytoplasmic proteins.
Proteome differences between treatments are linked with their functional roles
To understand the functional characteristics of the DEPs, Gene Ontology (GO) annotation was performed to analyze the genes and pathways involved in protein transformation. The GO terms of all DEPs were highly enriched in processes involved in catalytic activity, metabolic processes and cellular processes (Fig. 6A). Compared to those in the FR silage, the downregulated DEPs in the CK silage were enriched mainly in cellular process, metabolic process, and biosynthetic process; in contrast, the upregulated DEPs were enriched in catalytic activity, hydrolase activity, and peptidase activity (Fig. 6B, C). The upregulated processes related to protein transformation are shown in Figure 6D, including proteolysis, serine-type peptidase activity, endopeptidase activity, aspartic-type peptidase activity, etc. Figure 6E shows the presence of proteolytic enzymes in the silages, such as A0A396GND9 (putative actinidain), G7KU92 (putative tripeptidyl-peptidase II), G7I467 (putative nepenthesin), A0A396JDJ0 (carboxypeptidase), Q2HTQ3 (peptidase C1A, papain), G7IQG4 (fungal proteinase A, aspartic proteinase superfamily protein), G7KGD8 (putative prolyl oligopeptidase), A0A072U170 (subtilisin-like serine protease), etc. In the LN vs. CK group (Fig. 6F) and LB vs. CK group (Fig. 6H), most GO terms were downregulated, such as catalytic activity, metabolic process, hydrolase activity, etc. The processes related to protein transformation were also downregulated, e.g., peptidase activity, serine-type peptidase activity, and endopeptidase activity (Fig. 6G, I). These process variations were associated with downregulated levels of proteases, e.g., A0A396GND9 (putative actinidain), G7KGD8 (putative prolyl oligopeptidase), G7LCZ9 (serine carboxypeptidase S28 family protein), and G7J0H7 (aspartic proteinase nepenthesin-like protein) (Fig. 6F, H). Next, we investigated unique proteases between the LN and LB treatments. G7I467 (putative nepenthesin), G7J6G3 (Papain family cysteine protease), A0A072UV98 (Serine carboxypeptidase), and A0A072UWF5 (Eukaryotic aspartyl protease family protein) were uniquely detected in the LN vs. CK group, while G7ITL1 (Carboxypeptidase), G7KEU6 (Subtilisin-like serine protease), G7KU92 (putative tripeptidyl-peptidase II), and B7FKI4 (Carboxypeptidase) were exclusively detected in the LB vs. CK group (Fig. 6G, I). Functional analysis revealed that the ensiling process mainly altered proteins involved in cellular processes, catalytic activity, metabolic processes, and hydrolase activity, whereas the LN and LB treatments further downregulated these pathways.
The peptide and amino acid levels of alfalfa silage reveal protein transformation signatures
To identify the changes in protein hydrolysate during ensiling, the peptide abundance and composition were analyzed. In alfalfa silages, a total of 58 peptides were identified, comprising 9 tripeptides and 49 dipeptides (Fig. 7A). Compared with the FR and CK silages (Fig. 7A), the LN- and LB-treated silages exhibited an overall increase in the abundance of most peptides, except for Gamma-Glu-Leu and Gamma-Glu-Phe. The top 20 peptides with high abundance are shown in Figure 7B, among which the dipeptide Ile-Pro had the highest abundance, followed by Pro-Val, Val-Ile, and Val-Phe. To determine whether the endogenous peptides in silages have health-promoting value, PeptideRanker was used to predict the bioactivity of the peptides. A total of 25 peptides were labeled as bioactive (score > 0.5, Table S2). To further understand the underlying health benefits of the peptides, BIOPEP-UWM was used to determine their biological activity. Most of the dipeptides possess bioactivities, such as dipeptidyl peptidase (DPP) IV inhibitory activity and angiotensin-converting enzyme (ACE) inhibitory activity (Fig. 7C, Table S2). Some dipeptides, including Leu-Trp, His-Leu, and Ile-Arg, have antioxidative functions. Comparative analysis of peptide profiles revealed that LN- and LB-treated silages were enriched in bioactive peptides exhibiting DPP-IV inhibitory (DPP-IVi), ACE inhibitory, and antioxidant activities. This enrichment suggests a potential enhancement of health-promoting properties in silages subjected to LN and LB treatments.
Next, the variation in AA levels was assessed (Table S3). Compared to those in FR, the majority (80%) of the AAs in the LN- and LB-treated silages increased, except for Lys, Gln, His and Asp (Fig. 7D). However, compared with those in CK silages, the levels of Gln, Lys, and Asp in LB- and LN-treated silages significantly decreased, while the levels of L-Glu, Trp and Arg significantly increased (Fig. 7E). The results indicate that the levels of most AAs increased after ensiling, and the application of LB and LN further modulated specific AAs, particularly those belonging to the acidic and basic categories.
Microbial community alteration and microbial metabolism contribute to peptide levels
To explore whether the inoculation of L. plantarum mediates protein alterations during ensiling, AQ-16S-seq was conducted to assess microbial community changes. PCoA of the microbial community structure revealed significant differences between the groups (Fig. 8A). Figure 8B shows that the ensiling process and L. plantarum treatment resulted in a decrease in the community richness (Chao1, ACE, Shannon, and observed species), among which the LB treatment silages had the lowest richness. The relative abundance and absolute copy numbers of the identified microbiota at the species level are shown in Figure 8C. After ensiling, the absolute abundance of the total bacteria increased in the CK, LB and LN groups compared to that in the FR group. Compared with those in the CK group, both the relative and absolute abundances of Enterococcus mundtii and Bacillus sp. decreased, whereas the abundance of L. plantarum concurrently increased. The absolute abundance of L. plantarum was greatest in the LB group (9.11×107 copies/g sample, accounting for 98.76%), followed by the LN group (7.78×107 copies/g sample, accounting for 96.81%) (Fig. 8C). Linear discriminant analysis effect size (LDA effect size) revealed a distinct species, with L. plantarum having the highest LDA score, signifying its significant influence on the microbial composition (Fig. 8D, E). These results suggest that L. plantarum may be the key beneficial microbe that competes with undesirable microorganisms.
To reveal protein alteration mechanisms under L. plantarum treatment, functional protease genes were predicted according to KEGG databases. Through STAMP analysis, all identified protease genes were visualized with a heatmap (Fig. 8F), and the differences between groups were visualized with an extended error bar (Fig. 8G, H). Compared with those in the CK group, the levels of most protease genes were significantly lower under L. plantarum treatment, except for those of D-alanyl-D-alanine dipeptidase (in the LN and LB groups) and Pip (in the LN group). These results indicate that microbial enzymes may contribute to protein alterations during ensiling, whereas PrtP has a minimal impact.
Identification of the relationships between microbes and metabolites
To characterize the impact of functional species on the fermentation profiles, nitrogen availability, and protein fractions, Spearman’s rank correlation analysis was performed between distinct species and the aforementioned parameters. Strong positive correlations (correlation coefficient > 0.8, p < 0.05) between L. plantarum and LA/AA, Val-Ile-Arg, Trp-Phe, Ala-Val, Gly-Met, Leu-Pro-Ile, Ile-Leu, Pro-Val, Ile-Leu-Lys, His-Pro, and Val-Ile-Lys were detected. Conversely, there was a strong negative correlation (correlation coefficient < -0.8, p < 0.05) with AA-N, NH3-N, DM loss, Thr, His, Gln, acylaminoacyl-peptidase [EC:3.4.19.1], protease III [EC:3.4.24.55], pyroglutamyl-peptidase [EC:3.4.19.3], oligopeptidase A [EC:3.4.24.70], oligopeptidase B [EC:3.4.21.83], beta-aspartyl-dipeptidase (metallo-type) [EC:3.4.19.-], muramoyltetrapeptide carboxypeptidase [EC:3.4.17.13], and carboxypeptidase Taq [EC:3.4.17.19] (Fig. 9A, B). The correlation analysis results for E. mundtii, Bacillus sp., and W. paramesenteroides were opposite to those for L. plantarum. Furthermore, L. plantarum was also positively correlated with DM, PB2, peptide-N, Ile-Ile-Arg, etc., but inversely correlated with pH, ADF, NDF, etc. (Fig. 9C-G). In summary, L. plantarum contributed to a decrease in pH, LA formation, and the preservation of true protein and peptide-N, while E. mundtii, Bacillus sp., and W. paramesenteroides were associated with DM loss, NH3-N formation, and the waste of protein and peptides.