Short Communication: Prognostic Values of a Multiparametric Risk Score in Maize Silage Undergoing Different Ensiling Conditions

Abstract

We studied the effects of the use of Lactobacillus buchneri (Lb) and the maize pre-ensiling composition on the aerobic silage stability in a panel of 88 maize ensiled 60 days in 21 L buckets. Lb was dispensed at three dosages and compared to a control (pure water). The prognostic multiparametric risk score was used to find the risk factors related to the chemical composition of the fresh plant, associated with the onset of aerobic instability in maize silage. A multivariable Akaike’s Information Criterion in the backward Cox proportional hazard regression was estimated for pre-ensiled maize chemical traits. A Multiple Factorial Analysis (MFA) was calculated. The hazard ratios were 1.02, 1.34, 0.66, 0.65, 1.57, and 1.06 for dry matter (DM), crude protein (CP), ether extract (EE), aNDF, lignin (sulfuric acid, sa), and water-soluble carbohydrates (WSC), respectively (p < 0.05, DM, p = 0.15). At the MFA, ash, CP, aNDF, ADF, and lignin (sa) were grouped with a positive Dim-1, while DM, EE, and starch were grouped with a negative coordinate; WSC stood alone with Dim-1 close to zero. CP, EE, aNDF, lignin (sa), and WSC resulted in the most relevant traits and were used to build the nomogram. The use of strains of Lb improved the aerobic stability for maize harvested at <300 g/kg of DM.

1. Introduction

Maize silage represents the primary fiber and energy source in dairy cattle, but the aerobic deterioration, occurring typically in the feed-out phase [1], reduces its efficiency [2]. Maize is a valuable source of energy and nutrients and is feasibly incorporated in the total mixed ratio (TMR) [3]. The aerobic stability of maize silage is a consequence of several factors, such as the fermentative profile [4], the presence of oxygen during conservation [5], the ensiling practices [6], the environmental temperature during conservation [2], and the use of inoculants [7]. Silage deterioration can be favored by fungi, yeasts, and bacteria (primarily Bacillus species), consuming rapidly degradable carbohydrates, organic acids, and proteins, and producing heat [8,9,10]. Advanced technology permits the collection of big data from spatial yield and nutrient characteristics at harvest. These data can be used by farmers and maize growers to make better processing choices at the time of harvest and when processing the whole maize plant into the silos [11,12]. Due to big data availability, the use of Machine Learning (ML) techniques has become more common in agricultural research in establishing knowledge-based farming systems [13,14] and particularly for estimating silage fermentation quality from freshly harvested maize [15]. Maize is the most studied crop using ML because of its versatile usage and wide world diffusion, and even for its quality evaluation via the use of spectrometers [14]. Specifically, Near-Infrared Spectroscopy (NIR) instruments are known [16,17] and commercially available in harvesting machines.

The strains of Lactobacillus buchneri (Lb) increases acetic acid production in silage, inhibiting the yeasts and increasing aerobic stability [18], but the effects of different dosages on the fermentative profile are controversial [7,11,19].

The prognostic multiparametric risk score (PMRS) evidences the risk factors linked to the onset of a given disease, and nomograms are widely used because of their ability to reduce predictive models into a single numerical estimate of the probability of an event [20]. The PMRS was proposed to quickly evaluate the aerobic instability of maize silage, which can be considered a silage disease [11], and to find the risk factors associated with the onset of instability.

The objective of this study was to evaluate the risk of aerobic instability onset associated with some factors, such as the use of different dosages of Lb, the progressive delays in silos sealing, and the nutrient composition at harvest. All the studied factors can be monitored during harvesting by embedded sensors or controlled in the ensiling time. The chosen factors are more likely to be used to monitor and correct the silage production chain in agreement with the precision agriculture principles. Moreover, a nomogram to estimate the risk of silage degradation onset easily and rapidly was drawn.

2. Materials and Methods

We use Lb (CCM 1819) 4.2 × 10¹⁰ CFU/g (KWS LACTOSTABILITY, AGRAVIS Raiffeisen AG, Munster, Germany) in a panel of 88 maize samples of a late-ripening class hybrid (KWS KELINDOS, FAO class 600–130 d), sown in Northern Italy over three years (2018, n = 16; 2019, n = 36; 2020, n = 36), and harvested with a kernel unit and at a theoretical length cut of 12.7 mm. Lb was dispensed at three dosages (standard dose, SD = 2.02 × 10⁵ CFU/g of fresh forage, n = 26; half dose, HD = 1.01 × 10⁵ CFU/g of fresh forage, n = 18; double dose, DD = 4.04 × 10⁵ CFU/g of fresh forage, n = 18), and compared to a control (pure water, C, n = 26). Freshly harvested maize (FHM) was ensiled after 0 (D0, n = 40), 6 (D0, n = 24), or 20 (D0, n = 24) hours of air exposition. The 2018 samples only underwent the SD (n = 8), or C (n = 8), and D0. Samples were ensiled in a 21 L truncated conic plastic bucket and pressed using a 1-ton hydraulic press (141 kg/cm²) for 60 days. The buckets were shielded by using a 150 μm SealPlus Film permeable to oxygen at the rate of 48 cm³/m²/24 h at 23 °C and 65% RH (SealPlus by Gamma Srl, Mondovi CN Italy) and sealed with robust tape. The sealed buckets were stored in a dark room for 60 days at a stable temperature of approximately 23 °C, assuring a better anaerobic ensiling condition. The FHM was immediately analyzed for proximate composition using a portable NIR system (PoliSPECNIR, ITPhotonics, Breganze, Italy) [11,21]. The estimated chemical traits were: dry matter (DM), ash, crude protein (CP), fat (EE), aNDF, ADF, lignin (sulfuric acid, sa), water-soluble carbohydrates (WSC), and starch. Data are expressed as g/kg of the DM. Density was calculated as the ratio of the amount of pre-ensiled maize, corrected for the DM content, to the volume of the 21 L bucket, and was expressed as kg/m³. The porosity was calculated according to the formula proposed by Richard et al. (2004) [22].

With the intent of simulating the feed-out rate in a farm-scale silo-bunker, at opening time, at 48 and 96 h after silage opening, a 15 cm thick layer of silage was removed from the buckets, and a data logger was repositioned 7 cm under the silage surface, recording the temperature every 15 min with a precision of 0.1 °C (Elitech USB Temperature Datalogger RC-5, London, UK). The silage surface was cut by 10 cm and gently lifted on one side using a spatula; after introducing the data logger, the slice was repositioned, making sure to close the cut well. The event of aerobic instability for silages was defined as when their temperature exceeds 2 °C above room temperature [18,23]. Survival analysis, or more generally, time-to-event analysis, analyze the length of time until a well-defined endpoint, the event of interest, occurs [24]. For the survival analysis, the event was defined as the aerobic instability occurrence and the outcome as the time, expressed in hours, to the first event of aerobic instability onset (aerobic deterioration), and observations were censored at 108 h. A single data logger was used to record the room temperature throughout the entire period.

The normal distribution for continuous variables was assessed using the Shapiro-Wilk test and visual inspection. Preliminary descriptive statistics (Minimum, Maximum, Quartiles, Mean, and Standard Deviation values) were compiled for all variables of pre- ensiled datasets. A multivariable Akaike’s information criterion in the backward logistic regression (Log-AIC) was assessed on the predictive traits to estimate the event. Receiver Operating Characteristic (ROC) and Area Under the Curve (AUC) were used to check the classification model performances, and the confidence intervals (CI) were reported. The Survival curves (Kaplan-Meier, KM) were drawn, and log-rank calculated. A univariable Cox proportional hazard model (Cox-univariable) was estimated for the qualitative variables (the use of Lb, the sealing delays, and the years of harvesting). A multivariable Akaike’s Information Criterion in the backward Cox proportional hazard regression (Cox-AIC) was estimated for the pre-ensiled maize chemical traits. The proportional hazard assumption was evaluated using a visual approach (Schoenfeld, Martingala, beta, and score residuals) and Grambsh and Therneau test [25,26]. Drawing a nomogram allows a straightforward computation of the probability of occurrence of the event [20,27]. A Multiple Factorial Analysis (MFA) [28] was used to explore the relationship between pre-ensiled maize chemical traits, delay, and use of inoculant with the aerobic instability. The MFA is a multivariate data analysis used to sum up and simplify the data by reducing the dimensionality of the data set. The MFA considers different types of variables, types of structure of the data, and supplementary information. The use of R software version 4.0.2, 22 June 2020, performed all of the statistics, besides Rcmdr package version 2.6, and RStudio version 1.2.1578. The statistical significance was set at a p-value lower than 0.05.

3. Results

The pre-ensiled maize composition is reported in Table S1 of the Supplementary Materials. The effects of the years (2018, 2019, and 2020) or the outcome (no events = 0, events = 1), for pre-ensiled maize composition traits are reported in Table S2 of the Supplementary Materials. The number of samples recorded for the event was 31 of the initial 88. The Log-AIC showed AUC = 0.78 (95% C.I. = 0.67−0.89), AIC = 97.5, odds = 1.70, 0.64, 0.93, 1.37, 1.08 for the CP, EE, aNDF, lignin (sa), and WSC, respectively (p < 0.05, except for lignin (sa) with p = 0.12), Concordance index = 0.73.

Before drawing the KM curves, the sample silages were classified according to the DM content as follows: very low ≤ 300, 301 < low ≤ 330, 331 < medium ≤ 360, 361 < high ≤ 400, and very high > 400. The KM showed a tendential lower risk for maize inoculated with Lb when the DM at harvest was very low (p = 0.079) (Figure 1A). The KM showed differences among the years for inoculated (p = 0.05) (Figure 1B) or tendential differences among the years for non-inoculated samples (p = 0.10) (Figure 1C). Moreover, KM shows a tendential lower risk for inoculation with Lb harvested in 2019 (Figure 1D). The KM for the different years (2018, 2019, and 2020) stratified per Lb dosages, and the KM for the use of Lb, stratified per delays, are reported in Figures S1 and S2, respectively.

The Cox-univariate results are reported in Table S3 of the Supplementary Materials. The Cox-AIC, the Somer’s Dxy concordance index was 0.73, the Likelihood ratio test had p = 0.01, and the Schoenfeld residuals were graphically evaluated to ensure the proportional hazard assumption (global test had p = 0.17). The hazard ratios (HR) were 1.02, 1.34, 0.66, 0.65, 1.57, and 1.06 for DM, CP, EE, aNDF, lignin (sa), and WSC, respectively (p < 0.05, except for DM with p = 0.15). At the MFA, the first two latent variables, Dim-1 and Dim-2, explained 46.8 and 14.3% of the variability, respectively. DM, ash, EE, aNDF, ADF, lignin (sa), and starch had cos² > 0.5 in Dim-1, while WSC showed cos² > 0.5 in Dim-2. Ash, CP, aNDF, ADF, and lignin (sa) were grouped with a positive Dim-1, while DM, EE, and starch were grouped with a negative coordinate; WSC stood alone with Dim-1 close to zero. Along Dim-2, all traits showed positive values or close to zero, with the exception of WSC, showing a value of −0.83 (Figure 2A,B). At the nomogram (Figure 3A–C), the 95% C.I. risk for aerobic instability corresponded to a point range of 239–168 and 230–159 for 48 and 108 h, respectively.

4. Discussion

The occurrence of a well-defined event (i.e., the aerobic instability onset) is a binary outcome (0 or 1, stability or instability). This problem can be solved by Log regression, but if the time to the event is involved in the question, the survival analysis becomes more appropriate. The KM curves can graphically represent the plotting survival over time or estimate survival probability beyond a prespecified time interval. The Cox models relate the survival time to the covariate. Log and Cox models can be eventually adjusted (stratifying) for potential confounders. Those who survive until the end of the study or are lost in the follow-up have an unknown survival; instead, they are referred to as right censored [24]. The odds in Log and the HR in Cox models are the parameters that estimate the probability that the event occurs when comparing two groups.

The aerobic stability of maize silage depends on several factors, such as the use of Lb, the presence of oxygen, and a sufficient number of aerobic organisms (yeasts, molds, bacteria) and available substrates (sugars, starches, organic acids), a temperature > 4 °C and a DM content < 800 g/kg [29]. The obligate heterofermentative inoculants use lactic acid to produce acetic acid [7], primarily inhibiting the yeast activity responsible for the beginning of the aerobic spoilage [30]. The choice of the use of Lb among other possible species was related to its ability to increase the number of hours that a TMR remains stable when exposed to air and without any effects on DM intake or milk production [31]. In particular, one Lb among those commercially available was chosen. The higher DM content of the pre-ensiled maize could be a slight risk factor for aerobic instability [11], and this could be partially due to the lower WSC content available at later harvested maize, usable by the bacteria producing lactic and acetic acid [32], therefore, pH typically increases [29]. Here, at the Cox-AIC, the DM showed an HR of 1.04, confirming the finding of the previous study conducted by the authors [11]. The HR with values from 0 to 1 are preservative to risk, if equal to 1 expresses no effect, while if greater than 1 increases risk. The reported results in the Cox-AIC models show small-medium HR magnitudes and used factors that could vary when different algorithms for variable selection are used. The HRs for individual variables are not easily interpretable, but all together, they are representative of the plant maturity at harvest. However, we reported a difference in DM at harvest within the three experimental years. Samples inoculated with Lb showed a higher survival rate in 2018 and 2020 than in 2019 (Figure 1B), while in 2019, the inoculated samples showed a higher survival rate than the non-inoculated ones (Figure 1C). Although these results are only tendential, they are compatible with a previous study reporting the improved aerobic stability of silages treated with Lb compared to untreated ones, with a greater effect for samples with moderately high (410 g/kg) compared to normal (330 g/kg) DM at harvest [33]. However, within the three years studied, many unrecorded effects could have occurred due to the seasonality of climate; therefore, the years could be a confounder, and KM could benefit from the stratification by year. These findings suggest that the use of Lb improves aerobic stability in late harvest maize, and it plays a role even in wet samples.

Nevertheless, the different behaviors across the three years could not be explained by the average DM content alone. Indeed, conversely to the previous study [11], both Log-AIC and Cox-AIC indicate CP, lignin (sa), and WSC content as risk factors, while EE and aNDF are protective. These findings confirm the uncertain effect of the maturity stage of the plant at harvest on the aerobic instability onset [32], conversely to its well-known effect on the intensive rate of fermentation and acetic acid production during the ensiling process, and on the yeast counts content [34]. It can be speculated that other unknown variables may have occurred, e.g., the seasonality, within the aerobic stability process [32]. Moreover, the pre-ensiled chemical traits are reported to be significantly different in the three experimental years (Table S2). However, although the WSC content was similar among the three years, the aerobic instability samples showed higher WSC content.

Conversely, the other pre-ensiled chemical parameters showed differences among the years, but they had no differences in the averaged values between stable and unstable samples. Between the quantitative variables, neither density nor porosity were statistically associated with the event onset risk results, despite their variability among the test years (Table S2). However, due to the higher DM value in 2020, the density and porosity were higher (Table S2), probably due to a difficulty in compacting the maize mass into the silos. Furthermore, the porosity result was higher in 2020, proving a more consistent presence of air in the silage mass. However, recorded values for density (Table S1) were consistently lower than the recommended threshold of 240 kg/m³ or 700 as fed kg/m³, to limit the porosity at 0.40 [29]. These findings are confirmed by Table suppl mat 3, where it can be noticed that 2019 and 2020 were more at risk than 2018.

Despite the sealing delay being proven to influence the effectiveness of the inoculant treatment, to increase silage pH, and to affect the acetic acid [11,35,36,37], we did not find a role in the event’s occurrence (Figure S2, and Table S3). As for the use of Lb inoculant, which produced uncertain effects on aerobic stability (Figure 1, and Table S3), the difference in the dosages was unclear, reporting differences among the years for HD (Figure S2). This lack of results might be related to the deficiency of statistical differences from the use of HD, SD, and DD doses, while a log-scale difference in dosage would ensure appreciable differences [11].

The ML approach is advantageous in analyzing “big data” composed of highly variable and sample numbers, which are, unfortunately, limited in this study. Nevertheless, we wanted to explore the relationship between qualitative and quantitative predictors and outcomes. At the MFA (Figure 2A), there was no segregation between the samples that performed or not performed the event, and the qualitative variables had no trend among the two factors. Consistently with the MFA, the quantitative variables grouped for DM, EE and starch, and in the opposite, grouped for ash, CP, ADF, and lignin (sa), while WSC stands alone (Figure 2B). This finding confirms the positive relation between DM and starch [32], and data from the previous study, except WSC, where the spatial distribution of quantitative variables was studied with a Principal Component Analysis [15]. In the current study, we found an anomalous trend for the ratio starch/WSC during the three experimental years. However, the WSC can be related to both the maturity stage of the plant and some external factors such as cool, cloudy weather at harvest, or rainfall during wilting [29].

The use of pre-ensiled traits for PMRS could be implemented together with an NIR instrument directly in the harvesting machine. Nevertheless, the Cox-AIC model showed a concordance index of 0.73; therefore, more variables influence the event of aerobic instability, even in the case of a controlled ensiled process, such as in this study. Unfortunately, some of these variables cannot be monitored by NIR and, for such reason, were not included in the Cox-AIC predictive model used to build the nomogram.

However, despite the intention to provide a different approach in the study for the aerobic stability of maize, a more comprehensive dataset, including several variables, should be used to appreciate the potential of ML better, to explain the relationship between predictors and outcomes, and give a more stable definition of the PMRS approach.

As a limitation of the present study, it must be underlined that the results are obtained under experimental conditions and are applicable only to the studied factors. Therefore, validation using on-farm circumstances and variability is required to extend the results to the studied maize silage population.

Moreover, evaluation of the fermentative profile and the microbiological data was not in the scope of this trial but can improve the comprehension of the results and is recommended for future research. Although Lb is often used to strengthen silage’s aerobic stability, the use of different strains of Lb could lead to diverse aerobic stability results [38], and our findings are limited to the commercial Lb inoculated.

5. Conclusions

Pre-ensiled maize traits could be used as risk factors to predict the attitude of aerobic instability occurrence. Among the traits, CP, EE, aNDF, lignin, and WSC were the most relevant and were used to build the nomogram. Together with the use of a rapid sensor for maize characterization, these findings can be used in an automated system that will help farmers choose better strategies to enhance the quality of silage. The use of strains of Lb improved the aerobic stability for maize harvested at <300 g/kg of DM.