The fresh hybrid Pennisetum was harvested from a commercial plantation base located in Meizhou City, Guangdong Province (115.82°E, 24.52°N). Plants were cut down at a height of 2–2.5 m and chopped into pieces (1–2 cm per segment) via a fully automatic grass shredding machine. The chemical composition of hybrid Pennisetum was determined (Table 1). After moderate drying, the chopped hybrid Pennisetum was divided into five treatments and ensiled with (1) no additive (CK), (2) 1% malic acid (purity ≥ 99.5%; Shanghai Macklin Biochemical Co., Ltd., Shanghai) of fresh matter (MA), (3) 1% glucose (purity ≥ 99.5%; Shanghai Macklin Biochemical) of fresh matter (GL), (4) 100 U/g cellulase (VTR Bio-Tech Co., Ltd., Zhuhai, Guangdong) of fresh matter (CE), and (5) 10⁶ cfu/g Bacillus subtilis (VTR Bio-Tech) of fresh matter (BS). After adequate mixing, approximately 200 g hybrid Pennisetum was packed into a polyethylene bag (20 cm × 30 cm), compacted, and sealed via an automatic vacuum packager. Six replicates were set in a treatment; 30 bags were obtained and preserved at ambient temperature (25–28°C). After a 120-day fermentation, all 30 packages were opened for subsequent experiments, including the determination of fermentation parameters and chemical compositions and the analyses of bacterial community and metabolic profile.

Approximately 4 g of each silage sample was divided into two parts and maintained at −80°C for the determination of bacterial community and metabolic profile. Then, 10 g sample was homogenized with 90 ml distilled water and incubated at 4°C for 24 h. After extraction, the mixture was filtered by a sterilized 4-layer gauze, and the filtrate was collected to measure fermentation parameters. pH was determined using a glass electrode pH meter (FE28-Standard); volatile fatty acids (VFAs) and LA were determined in accordance with the method of Rumsey et al. (1967); and ammoniacal nitrogen (AN) was measured via phenol–hypochlorite colorimetry (Zhang et al., 2021).

The rest of the silage samples was dried to constant weight at 65°C for 48 h for the determination of dry matter (DM). Then, the dried sample was ground for the analyses of chemical compositions. WSC and crude protein (CP) contents were determined in accordance with the procedure of the Association of Official Analytical Chemists (AOAC, 2002). Neutral (NDF) and acid (ADF) detergent fibers were measured in accordance with the methods described by Van Soest et al. (1991). True protein (TP) and non-protein nitrogen (NPN) contents were determined using the trichloroacetic acid method (Licitra et al., 1996). Neutral (NDIN) and acid (ADIN) detergent insoluble nitrogen contents were measured in accordance with the method of Licitra et al. (1996). The content of starch was analyzed via perchloric acid–anthrone colorimetry reported by Bakhshy et al. (2020).

The total genomic DNA of silage sample was isolated via the DNeasy Power Soil Kit (QIAGEN, Inc., Netherlands) on the basis of the manufacturer’s instructions. After isolation, the concentration and quality of DNA samples were determined using the NanoDrop2000 Spectrophotometer (Thermo, United States). The Pyrobest DNA Polymerase (TaKaRa, DR500A) was adopted for the amplification of 16S rRNA V3–V4 regions of genomic DNA, and primer pairs were designed as 338F (5′-ACTCCTACGGGAGGCA GCA-3′) and 806R (5′-GACTACHVGGGTATCTAATCC-3′) in accordance with the method of Zi et al. (2021). The quality control and purification of PCR products were subsequently performed after amplification via the PicoGreen dsDNA Assay Kit (Invitrogen, Carlsbad, CA, United States) and Agencourt AMPure Beads (Beckman Coulter, Indianapolis, IN, United States). The high-throughput sequencing was carried out by the equimolar and paired-end sequencing (PE250) on the Illumina Novaseq 6000 platform (Personal Biotechnology Co., Ltd., Shanghai, China).

After sequencing, the processing of sequenced data was implemented by the QIIME (V 1.8.0) software. After the filtration of low-quality sequences defined by Gill et al. (2006) and Chen and Jiang (2014) and meaningless sequences (Including adapters, chimera, poly_A and primer), valid sequences were clustered into operational taxonomic units (OTUs) at a 97% similarity threshold via the UCLUST (Malik et al., 2020). The Basic Local Alignment Search Tool was used for further taxonomic classification, and an OTU table was obtained (Urbanek et al., 2020). The α- and β-diversity values of bacterial community were calculated via the QIIME software and vegan package in R software, respectively (De Filippis et al., 2018). The functions of bacterial community were predicted via the Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) database (Langille et al., 2013). The Linear Discriminant Analysis (LDA) Effect Size (LEfSe) analyses were conducted via an online tool¹, and LDA score>3 and p < 0.05 were selected as threshold.

In this study, the metabolic profile was analyzed via the UHPLC-Q/TOF-MS technology. The thawed sample (100 mg) was ground promptly and homogenized in 1 ml precooled methanol/acetonitrile/ddH₂O solvent (2:2:1, v/v/v), and the mixture was subjected to cryogenic ultrasound at −20°C for 30 min. After standing for 10 min, the mixture was centrifuged at 14,000 rpm and 4°C for 15 min. The supernatant was reserved and dried by a vacuum centrifuge, and the dried sample was redissolved and homogenized in 100 μl acetonitrile/ddH₂O solvent (1:1, v/v). After centrifugation at 14,000 rpm and 4°C for 15 min, the supernatant was collected for analyses. In addition, quality control (QC) samples were prepared to monitor the repeatability and stability of instruments.

Samples were analyzed using the Agilent 1290 Infinity LC UHPLC system and the parameters of instrument were set as follows: column temperature, 25°C; flow velocity, 0.5 ml/min; injection volume, 2 μl; mobile phases, and 25 mM ammonium acetate + 25 mM ammonium hydroxide in ddH₂O (A) and acetonitrile (B). The elution procedure was as follows: 0–0.5 min, B maintained at 95%; 0.5–7 min, B changed from 95 to 65% linearly; 7–8 min, B changed from 65 to 40% linearly; 8–9 min, B was 40%; 9–9.1 min, B linearly changed from 40 to 95%; 9.1–12 min, B maintained at 95%. Samples were stored in an autosampler at 4°C during the whole analysis period.

The AB Sciex Triple TOF 6600 mass spectrometer was adopted for analyses, and the parameters of electrospray ionization source were set in accordance with the method of Li R. et al. (2021). Briefly, the Ion Source Gas1, Ion Source Gas2, and Curtain gas were 60, 60, and 30 kPa, respectively. The IonSpray Voltage Floating was ±5,500 V for positive and negative modes. The source temperature was 600°C. For MS-only acquisition, the m/z range of the TOF MS scan and product ion scan were 60–1,000 and 25–1,000 Da, respectively. The accumulation time of TOF MS scan and product ion scan were 0.20 and 0.05 s/spectra, respectively. The information-dependent acquisition was adopted for the acquisition of the product ion scan, and parameters were set as follows: mode, high sensitivity; declustering potential, ±60 V for positive and negative modes; collision energy, 35 V ± 15 eV. Ten candidate ions were supervised per cycle, and isotopes within 4 Da were eliminated.

Data processing was performed in accordance with the method of Li R. et al. (2021). First, raw data in the wiff.scan format were converted into the MzXML format via the ProteoWizard MSConvert, and MzXML files were imported into the XCMS software. The parameters for peak picking were set as centWave m/z = 25 ppm, prefilter = c (10, 100), and peak width = c (10, 60). bw = 5, minfrac = 0.5, and mzwid = 0.025 were used as parameters for peak grouping. The annotation of adducts and isotopes was performed using the Collection of Algorithms of MEtabolite pRofile Annotation. The structural identification of metabolites was performed by comparing the retention time, molecular weight (<25 ppm), MS/MS spectra, and collision energy with the database established by Benton et al. (2015). Subsequently, multidimensional statistical analyses, including principal component analysis (PCA), projections to latent structures–discriminant analysis (PLS-DA), and orthogonal PLS-DA (OPLS-DA), were conducted using the SIMCA-P 14.1 software. Metabolites with variable importance in projection (VIP) > 1, fold change > 1.2 or <0.833, and p < 0.05 were recognized as differential metabolites (Hao et al., 2021). In addition, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database was adopted for the analyses of functional characteristics and classifications of differential metabolites (Guan et al., 2020).

The experimental data of fermentation parameters, chemical compositions, and α-diversity indices were preliminarily sorted using the Excel software and analyzed via the SAS 9.4 software with the model of Y_ij = μ + A_i + ε_ij. Y_ij is the dependent variable of the silage samples in different treatments, μ is the overall mean, A_i is the effects of silage additives, and ε_ij is the random error. The LSD method was adopted for multiple comparisons. Data were displayed in tables in the form of mean ± SD, and p < 0.05 indicated a significant difference. The Pearson correlation coefficient was adopted for correlation analyses, and correlation coefficient (Cor) > 0.6 or <−0.6 and p < 0.05 were considered as correlation.

The chemical compositions of the hybrid Pennisetum are shown as Table 1. The DM content of hybrid Pennisetum was 336.30 g/kg fresh matter (FM). The contents of WSC, starch, NDF, ADF, and CP were 26.84, 8.78, 696.50, 431.72, and 68.05 g/kg DM, respectively. The contents of TP, NPN, NDIN, and ADIN were 742.91, 257.09, 252.61, and 188.91 g/kg total nitrogen (TN), respectively.

As shown in Table 2, glucose and cellulase increased and decreased, respectively, the DM content in HPS (p < 0.05). The application of malic acid and cellulase decreased the pH of HPS (p < 0.05). Malic acid increased the content of LA but decreased the AA content (p < 0.05). Compared with that in CK, the AA content in CE continued to increase (p < 0.05). WSC contents in GL, CE, and BS were significantly lower than those in CK (p < 0.05). The contents of AN in MA and CE decreased compared with that in CK (P < 0.05). The chemical compositions of HPS are shown in Table 3. Malic acid, cellulase, and B. subtilis reduced the contents of NDF in HPS (p < 0.05). Moreover, cellulase and B. subtilis could reduce the contents of starch and ADF (p < 0.05). Compared with CK, GL, CE, and BS were observed with increased CP and significantly decreased NDIN and ADIN (p < 0.05). In addition, the application of malic acid and cellulase could increase the content of TP (p < 0.05).

In this study, silage samples were sequenced via the Illumina Novaseq 6000 platform, and 2,548,477 raw reads were obtained. After removing meaningless sequences, we obtained 2,011,503 clean sequences. The clustering and annotation of clean reads were performed, and 13,503 OTUs were found eventually. The diversity of bacterial communities is shown in Figure 1. The rank–abundance curve showed that samples in MA had improved evenness and diversity (Figure 1A). With increasing sequence number, the rarefaction curve tended to flatten, indicating that sampling was sufficient and that data were reasonable and reliable for further analyses. For α-diversity (Figures 1C–H and Supplementary Table 1), the Observed species, Faith PD, Pielou evenness, Chao1, Shannon, and Simpson indices of MA were significantly higher than those in CK (p < 0.05). Moreover, the Pielou evenness and Simpson indices in CE were lower compared with those in CK (p < 0.05). The Venn diagram showed that 351 OTUs appeared in all five treatments, and 1,678, 3,389, 1,423, 2,259, and 1,969 OTUs were unique to CK, MA, GL, CE, and BS, respectively (Figure 1I). For β-diversity, the PCoA score plot indicated that the individuals in MA could be remarkably separated from those in four other treatments. The Adonis analyses were further conducted, and results showed that the application of malic acid, cellulase, and B. subtilis had significant effects on the bacterial communities in HPS (p < 0.05). However, no evident difference between CK and GL was observed (p < 0.05).

At the phylum level, the dominant bacteria were Proteobacteria, Firmicutes, Bacteroidetes, Actinobacteria, and Deinococcus–Thermus (Figure 2A and Supplementary Table 2). Malic acid reduced the relative abundance of Proteobacteria but increased Firmicutes and Actinobacteria (p < 0.05). The addition of B. subtilis increased the abundance of Deinococcus–Thermus (p < 0.05, Figure 2B and Supplementary Table 2). At the genus level, the dominant bacteria were Aquabacterium, Bacillus, Weissella, Lactobacillus, and Brevibacillus (Figure 2A and Supplementary Table 2). The application of malic acid decreased the abundance of Aquabacterium and Weissella but increased the abundance of Bacillus and Brevibacillus (p < 0.05, Figure 2B and Supplementary Table 2).

To study the differential OTUs and potential biomarkers, we performed LEfSe analyses between CK and each additive treatment. As shown in Figure 3 and Supplementary Table 3, the differential OTUs between CK and MA were highest, and the addition of malic acid caused the upregulation of 81 differential OTUs and the downregulation of 48 differential OTUs (Figure 3A). Glucose caused the least number of differential OTUs, and the numbers of upregulated and downregulated OTUs were 3 and 18, respectively (Figure 3C). The number of differential OTUs caused by cellulase was 8 for upregulation and 17 for downregulation, and 17 upregulated and 23 downregulated OTUs were caused by B. subtilis (Figures 3B,D).

In this study, the PICRUSt was adopted for the functional prediction of bacterial communities (Figure 4 and Supplementary Table 4). The top five predicted functions were DNA helicase, DNA-directed DNA polymerase, histidine kinase, NADH:ubiquinone reductase (H⁺-translocating), and peptidylprolyl isomerase (Figure 4A). The top five pathways were aerobic respiration I (cytochrome c), gondoate biosynthesis (anaerobic), CDP–diacylglycerol biosynthesis I, CDP–diacylglycerol biosynthesis II, and cis-vaccenate biosynthesis (Figure 4C). The diversity of functions and pathways were also calculated, and the PCoA score plots of functions and pathways indicated that treated individuals were separated from four other treatments. The Adonis analyses further confirmed that malic acid played important roles in altering microbial functions (Table 4), but the effects caused by glucose, cellulase, and B. subtilis were unnoticeable.

The top five differential pathways between CK and four additive treatments are shown in Figure 5 and Supplementary Table 5. For malic acid and B. subtilis, the application of the two additives reduced the abundance of the top five differential pathways (p < 0.05), whereas the abundance of pathways were upregulated by cellulase (p < 0.05). Glucose upregulated the abundance of the superpathway of thiamin diphosphate biosynthesis II but downregulated the abundance of the superpathway of pyridoxal 5′-phosphate biosynthesis and salvage, PpGpp biosynthesis, pyridoxal 5′-phosphate biosynthesis I, and Syringate degradation (p < 0.05). In addition, we could find that the abundance of differential pathways in CK vs. GL were less than those in CK vs. MA, CK vs. CE, and CK vs. BS, implying that the effects of glucose on the dominant pathways were not evident. By analyzing all differential pathways, we found that nine pathways were simultaneously affected by all four additives (Figure 6 and Supplementary Table 6).

The relative standard deviations (RSD) of the ion peak abundance in QC samples are displayed in Supplementary Figure 1A (positive ion mode) and Supplementary Figure 1B (negative ion mode). In this study, the number of peaks with RSD ≤ 30% in QC samples accounted for more than 80% of the total number of peaks, indicating that the instrumental analysis system was stable and reliable (Shen et al., 2016). The PCA score plots of 30 samples and 4 QC samples are shown in Supplementary Figure 2A (positive ion mode) and Supplementary Figure 2B (negative ion mode). Results showed that QC samples were closely clustered, suggesting that the repeatability of the experiment was good.

In this study, 422 metabolites, including 259 in the positive ion mode and 163 in the negative ion mode, were obtained. These metabolites could be divided into organic acids and derivatives, organic oxygen compounds, lipids and lipid-like molecules, organoheterocyclic compounds, benzenoids, and five other types (Supplementary Figure 3). Multivariate statistical analyses were conducted to distinguish the difference between CK and the four other treatments. PCA and PLS-DA score plots showed that the effects of malic acid and cellulase on the metabolic profile in HPS were evident, but the differences between CK and GL and between CK and BS were not significant (Figure 7). Subsequently, we performed OPLS-DA analyses, and the VIP values of each metabolite were calculated. After screening with the criterion of VIP > 1, FC > 1.2 or <0.833, and p < 0.05, differential metabolites were obtained (Supplementary Table 7). Malic acid upregulated 29 metabolites and downregulated 18 metabolites, and glucose upregulated 8 metabolites and downregulated 5 metabolites. The numbers of upregulated and downregulated metabolites caused by cellulase were 35 and 12, respectively. The numbers of upregulated and downregulated metabolites caused by B. subtilis were 16 and 2, respectively (Figure 7).

We further performed the functional annotation of differential metabolites via the KEGG database (Figure 8 and Supplementary Table 8). Results showed that differential metabolites were involved in carbon fixation in photosynthetic organisms; alanine, aspartate, and glutamate metabolism; aminoacyl-tRNA biosynthesis; glyoxylate and dicarboxylate metabolism; ascorbate and aldarate metabolism; beta-alanine metabolism; arginine biosynthesis; citrate cycle (TCA cycle); isoquinoline alkaloid biosynthesis; tyrosine metabolism; and betalain biosynthesis.

In the comparison between CK and MA, g_Thermobacillus was positively correlated with alanine–glucose (Cor = 0.9793, p = 2.92E−08) but negatively correlated with L-alanine (Cor = −0.9510, p = 2.05E−06) and D-threitol (Cor = −0.9419, p = 4.71E−06; Figure 9A). For CK vs. GL, s_Alistipes finegoldii was positively correlated with isomaltose (Cor = 0.8647, p = 5.92E−04), and f_Rikenellaceae was negatively correlated with alpha-D-glucose (Cor = −0.8351, p = 1.38E−03) and N-acetyl-L-phenylalanine (Cor = −0.7928, p = 3.60E−03; Figure 9B). In the comparison between CK and CE, s_Bacillus smithii was positively correlated with threonine–valine (Cor = 0.9085, p = 4.32E−05) and 5-methylcytosine (Cor = 0.8798, p = 1.61E−04), and s_Lactobacillus brevis was negatively correlated with isoleucine–isoleucine (Cor = −0.7959, p = 1.96E−03) and proline–glucose (Cor = −0.7741, p = 3.12E−03, Figure 9C). For CK vs. BS, s_Brevibacillus thermoruber was positively correlated with histidine–proline (Cor = 0.8935, p = 8.98E−05) and isomaltose (Cor = 0.8850, p = 1.30E−04), and tyramine was negatively correlated with g_Thermoactinomyces and f_Thermoactinomycetaceae (Cor = −0.9297, p = 1.20E−05; Cor = −0.9232, p = 1.85E−05; Figure 9D). Other correlations are displayed in Supplementary Table 9.

In this study, the DM content of the hybrid Pennisetum was 336.30 g/kg FM, which met the ideal DM content of 30–35% (Guyader et al., 2018). NDF, ADF, and CP contents were all higher than that reported by Shah et al. (2020a). The variation might be attributed to the geographical and climatic conditions for plant growth and the growth stage of a plant. The WSC content of hybrid Pennisetum was 26.84 g/kg DM, which was insufficient to meet the minimum requirement of 60–70 g/kg DM for well-preserved silage (Wang et al., 2021). As a result, the moisture content of hybrid Pennisetum was suitable for silage, but the insufficient WSC might prevent continuous fermentation.

Dry matter losses are caused by the metabolic activities of aerobic microorganisms, and digestible carbohydrates are an important fraction for consumption (Yuan et al., 2017). The addition of glucose might alleviate the consumption of carbohydrates to some extent. The pH values of HPS in all five treatments were less than 4.2 and could be considered as well fermented (He et al., 2020). Malic acid reduced the pH due to its acidity, and the decreased pH in CE might be attributed to increased AA. Cellulase can hydrolyze structural carbohydrates and provide a fermentation substrate for LAB, thereby accumulating organic acids and reducing pH (Li F. et al., 2020). The changes in the contents of LA and AA in CK and MA might be because of different fermentation types. WSC was consumed and decomposed into organic acids under the activities of microorganisms during the first 30 days of silage, and the fermentation in GL, CE, and BS might be complete (Chen et al., 2020). Similar to WSC, starch is a type of rapidly degradable carbohydrate and can provide a substrate for fermentation (Feng et al., 2020). The CP content in CE was highest, which might be due to the carbon loss during fermentation (He et al., 2018). In addition, silage additives had remarkable effects on the nitrogen fractions of HPS. AN is produced by clostridial fermentation, and one of the criteria for good silage is AN < 100 g/kg TN (Mu et al., 2020). NPN is an important product of proteolysis, and the utilization rate of NPN in ruminants is lower compared with that of TP due to rapid degradation (Dentinho et al., 2019; Wang et al., 2021). The present study suggested that the addition of malic acid and cellulase had positive effects on inhibiting proteolysis, and low pH was speculated to inhibit the activities of proteases (He et al., 2020). NDIN and ADIN belong to bonding proteins, which bond plants with cell wall. NDIN is degraded and utilized in the rumen slowly, whereas ADIN cannot be utilized (Licitra et al., 1996). As a result, the four additives used in this study had positive effects on the nutrients of HPS, and this finding was manifested in increased available protein and reduced structural carbohydrates.

Given the important influence of microorganisms on silage fermentation, bacterial communities should be monitored. The α-diversity refers to the richness and diversity of species in a particular habitat or ecosystem and can indicate the degree of species isolation in the habitat. According to Bai et al. (2020), when dominant bacteria are abundant, the diversity of microbial community is low. In this study, the dominant bacteria at the phylum (Proteobacteria) and genus (Aquabacterium) levels in CE were highest among the five treatments, whereas the dominant bacteria in MA were lowest. These findings might be the reason why the α-diversity indices in CK were lower and higher than those in MA and CE, respectively. The β-diversity can reflect the differences between and within groups. Results indicated that malic acid, cellulase, and B. subtilis had significant effects on the bacterial communities of the HPS. However, the effects of glucose seemed inconspicuous.

Bacterial communities were altered by silage additives, whereas Firmicutes and Proteobacteria were always the dominant phyla. Results were consistent with those reported by Wang et al. (2021); Zhang et al. (2021), and Zi et al. (2021). Almost all LAB belong to Firmicutes, which are involved in the degradation of biomacromolecules and secretes various lipases, cellulases, and proteases (Yuan et al., 2020). The upregulation of Firmicutes in MA might lead to increased LA. A wide variety of pathogenic Gram-negative bacteria, including Escherichia, Helicobacter, Legionellales, Salmonella, Vibrio, and Yersinia, belong to Proteobacteria. Malic acid may inhibit pathogenic bacteria (Du et al., 2021). Moreover, Actinobacteria and Deinococcus–Thermus were upregulated by malic acid and cellulase, respectively. However, research on the two phyla in silage remains limited, and mechanisms need further research.

At the genus level, Aquabacterium was the dominant bacteria and accounted for the vast majority of Proteobacteria. According to previous reports, Aquabacterium is a facultatively anaerobic, non-spore-forming, and rod-shaped bacteria detected in soil and fresh water, and most species cannot metabolize carbohydrates (Chen et al., 2012; Hirose et al., 2020). Thus far, only Xu et al. (2020) reported Aquabacterium in silage. Therefore, the function of the genus remains unknown. Bacillus, a type of facultatively anaerobic and Gram-positive bacteria, can secrete bacteriocin and inhibit the activities of undesirable bacteria. In addition, some species of Bacillus can produce LA (Li M. et al., 2020). Weissella, a type of heterofermentative LAB, converts 1 mol glucose to 1 mol AA, 1 mol LA, and 1 mol CO₂ (Wang M. et al., 2019). The upregulation of Bacillus and downregulation of Weissella may lead to increased LA and decreased AA in MA. Some species of Brevibacillus are reported to produce antibacterial agents (Miljkovic et al., 2019). The antimicrobial activities of malic acid might be attributed to the upregulation of Bacillus, Brevibacillus, and other similar bacteria.

The LEfSe analysis is a tool for the interpretation and discovery of high-dimensional data biomarkers, highlights biological correlation and statistical significance, and can explore potential biomarkers that are different among groups. In this study, only two differential bacteria were affected by all four additives. B. smithii was upregulated, whereas Lactobacillus rossiae was downregulated, indicating that the four additives had different regulating mechanisms on silage quality. B. smithii can convert cellulose into LA via simultaneous saccharification and fermentation processes (Chacón et al., 2021). Although research about B. smithii on silage remains limited, its characteristics of simultaneous saccharification and fermentation allow increased fermentation efficiency. B. smithii may be developed as a new silage or food additive, but the safety of produced silage or food should be assessed. L. rossiae, a type of LAB, converts arginine into toxic putrescine via the ornithine decarboxylase pathway (Del Rio et al., 2018). L. rossiae may be a biomarker of silage safety and quality, but no direct evidence shows that L. rossiae is related to silage spoilage and accumulation of toxins. Notably, Moraxellaceae is inhibited by glucose, cellulase, and B. subtilis and previously reported as a pathogenic microorganism (Narciso-da-Rocha et al., 2018). Lysinibacillus and Paenibacillus were upregulated by malic acid. These bacilli are reported to be probiotics and considered as potential fermentation additives (Liñan-Vidriales et al., 2021; Soni and Keharia, 2021). Methylobacterium is a strictly aerobic bacterium and involved in environmental carbon cycle and decomposition of compounds in plants (Wang et al., 2021). Ali et al. (2020) reported that Methylobacterium is downregulated from the pre-ensiled to the post-ensiled period, and the reduction of the genus implies that malic acid may create good conditions for fermentation. Leuconostoc is a lactate-producing bacteria (Wang et al., 2021). However, with increasing LA in MA, the abundance of the genus declined. This finding was speculated to be caused by competition between LAB, and the growth of dominant LAB might inhibit the metabolic activities of LAB with low abundance. The addition of glucose reduced the relative abundance of some pathogenic bacteria, including Escherichia–Shigella and Acinetobacter (Liu et al., 2020; De Almeida et al., 2021). Similarly, Acinetobacter was downregulated by B. subtilis. Pseudomonadales was downregulated by cellulase. According to a previous report, some members of the order are spoilage microorganisms with lipolytic and proteolytic enzymatic activities (Porcellato et al., 2018). This finding might explain why CE had the highest CP and TP contents among the five treatments. In summary, the results of LEfSe analyses implied that bacterial communities might be optimized by four additives. However, the mechanisms required in-depth research because of the complex competition and cooperation in microorganisms.

The function of bacterial communities was determined by their composition to a high extent. In this study, malic acid had significant effects on the functions of bacterial communities, as evidenced by the addition of malic acid caused the most differential bacteria compared with the addition of glucose, cellulase, and B. subtilis. This finding might be because the malic acid additive led to abundant variations in some functional bacteria (Wang et al., 2021). Malic acid has antioxidant and antibacterial functions and is widely used in food, pharmaceutical, healthcare, cosmetics, and other industries. Moreover, malic acid, which can be used as a carbohydrate source, provides energy for LAB and accelerates LAB growth (Li M. et al., 2020). Notably, the relative abundance of differential pathways caused by glucose was lower compared with those of pathways caused by malic acid, cellulase, and B. subtilis. Thus, glucose had little effect on dominant pathways. Nine pathways, including 1 upregulated (i.e., superpathway of thiamin diphosphate biosynthesis II) and 8 downregulated (i.e., PpGpp biosynthesis, syringate degradation, protocatechuate degradation I [meta-cleavage pathway], methyl gallate degradation, gallate degradation II, superpathway of salicylate degradation, toluene degradation IV [aerobic, via catechol], and nitrate reduction I [denitrification]) pathways, were affected by all four additives. We hypothesized that the nine pathways were the core pathways affecting silage quality in the hybrid Pennisetum.

Ensiling is a complex biological fermentation process that involves a wide variety of microorganisms, thereby leading to many kinds of metabolites. Conventionally, VFAs and LA are adopted for evaluating silage quality, whereas other metabolites are also of concern. Therefore, metabolomics technology can provide a complete understanding of metabolites in the microenvironment. The addition of glucose caused 13 differential metabolites, which was the least among the four additives. CK and GL could not be separated in PCA and PLS-DA score plots, which might be because glucose had less effect on the microflora and caused few differential metabolites. No metabolite was simultaneously upregulated or downregulated by all four additives. This finding was consistent with the fact that the four additives shared few differential bacteria and further suggested that the mechanisms of the four additives altering silage quality varied. In addition, the functional annotation based on KEGG database showed that the differential metabolites were involved in the metabolism and biosynthesis of small molecules, TCA cycle, and carbon fixation in photosynthetic organisms. This was beneficial to the comprehension of metabolites in silage and contributed to the identification of beneficial metabolites (Guan et al., 2020).

A total of 26 metabolites were obtained in MA vs. CK and CE vs. CK, indicating that malic acid and cellulase had the most similar effects on the metabolic profiles of HPS. Besides, malic acid and cellulase affected four pathways, including carbon fixation in photosynthetic organisms; alanine, aspartate, and glutamate metabolism; ascorbate and aldarate metabolism; and beta-alanine metabolism. Considering that malic acid and cellulase reduced the ratio of AN and NPN in HPS, we speculated that the four pathways might play some unknown functions. Among the four pathways, five differential metabolites, including L-alanine, L-aspartate, L-gulono-1,4-lactone, myo-inositol, and uracil, were obtained. In this study, L-aspartate was involved in five pathways and negatively related to B. smithii. Coincidentally, this species was negatively related to L-alanine, which played important roles in three pathways. This result confirmed the important role of B. smithii in silage. In the research of metabolome on silage, amino acids as metabolites are reported, whereas the functions of amino acids are still unknown (Guo et al., 2018; Xu et al., 2019). Besides, we found that malic acid and cellulase could upregulate the expression of several dipeptides, like alanine–glucose, threonine–glucose, histidine–isoleucine, and proline–glucose. Some dipeptides are reported to exert anti-inflammatory functions and considered as flavor substances (Gallego et al., 2021; Molinari et al., 2021). Thus, silage supplemented with malic acid and cellulase might provide improved palatability and animal health. Tuberibacillus, B. smithii, and Bacillus lehensis might play positive roles in the production of dipeptides. In addition, as an important part of the TCA cycle, malic acid can directly regulate the metabolism of carbohydrates and proteins through this biochemical reaction (Li M. et al., 2020). Two differential metabolites (i.e., malic acid and isocitrate) caused by glucose are involved in the TCA cycle and glyoxylate and dicarboxylate metabolism, which are reported to be related with carbohydrate metabolism (Ruan et al., 2021). According to the correlation analysis, bacteria associated with malic and isocitrate included Pseudomonadales, Acinetobacter, Pelomonas, and Lactobacillus reuteri, suggesting that these microorganisms might be potential targets for regulating carbohydrate metabolism during silage. For other metabolites, some special sugars are considered as biomarkers of poor fermentation. According to Li M. et al. (2021), the downregulation of isomaltose and L-sorbose may indicate high-quality ensiling. B. subtilis caused 18 differential metabolites. Among these metabolites, diosmetin is a flavonoid and has antioxidant, antibacterial, and anti-inflammatory properties (Lee et al., 2020). In the present study, diosmetin was positively related to Thermoactinomyces daqus and Brevibacillus thermoruber, suggesting that these two species might be involved in the production of diosmetin and other bioactive compounds. Several studies proved that different additives in silage can cause changes in bioactive metabolites, such as polyphenols, flavonoids, and terpenoids (Guan et al., 2020; Hu et al., 2020; Xu et al., 2021). However, these active substances are rarely used directly as silage additives. The effects of active substances on silage should be studied. Arachidonic acid (cis-5,8,11,14-eicosapentaenoic acid) is a type of ω-6 polyunsaturated fatty acid with anti-inflammatory functions (Trostchansky et al., 2021). Arachidonic acid was positively related to Bacillus, and the upregulation might be attributed to the inoculation of B. subtilis. Notably, correlation analyses were performed in accordance with statistical and correlation parameters. Therefore, results could not be regarded as causation (Xu et al., 2019). The present study provided potential biomarkers for improving the silage quality, and the development and application of amino acids, organic acid, small peptides, phenols, and other bioactive compounds will be the research focus of silage additives (Xu et al., 2019).