Using null models to decipher bacterial assembly mechanisms in oat silages harvested from southern China

Background Deciphering the assembly rules of microbial communities is vital for a mechanistic understanding of the general principles driving microbiome structures and functions. In this study, a null modeling-based framework was implemented to infer the assembly rules of bacterial community in oat silages harvested in southern China starting from the grain-filling stage through to full ripening. Results Most silages displayed “inferior” or “very inferior” fermentation quality. The fermentation qualities of silages tended to further decrease with the delay of harvest. Lactobacillus , Pediococcus , unclassified_f_Enterobacteriaceae , and Hafnia–Obesumbacterium constituted the predominated genera in silages. Delaying harvest increased the proportions of Hafnia–Obesumbacterium . Null model analysis revealed that stochastic processes were the primary contributor to the assembly of rare subcommunity during silage fermentation. The succession of abundant subcom-munity was controlled both by stochastic and deterministic processes. Deterministic processes, more specifically, heterogeneous selection, were more prominent in the assembly of abundant bacteria in silages with the delay of harvest. Linear regression analysis indicated the important roles of DM, WSC and pH in the assembly of abundant subcommunity. Conclusion This study, from the ecological perspectives, revealed the ecological processes controlling the bacterial community assembly in silage, providing new insights into the mechanisms underlying the construction of silage bacterial community.


Introduction
Forage production is seasonal in many parts of the world.To make the greens available throughout the year, ensiling is a common way to preserve the green fodder.Ensiling is an anaerobic bacterial-based fermentation process.During fermentation, sugars are converted by epiphytic lactic acid bacteria (LAB) mainly to lactic acid and acetic acid.This lowers the pH and creates an environment where the resulting silage is preserved.However, biochemical, and microbiological incidents can arise at the different stages of ensiling which result in high variability in fermentation quality.Poorly fermented silage would induce economic losses, affect animal performance, and even threaten animal and human health [1].
Silage fermentation process involves a variety of bacterial communities and biochemical reactions.Existing studies on deciphering silage fermentation mainly depends on the dynamics of bacterial community [2,3].Nevertheless, what ecological mechanisms govern the formation and development of bacterial community structures in silage is poorly known.From the ecological perspectives, microorganisms establish communities according to deterministic or stochastic processes.Deterministic theories suggest that local, niche-based processes, such as environmental filtering, biotic interactions and interspecific trade-offs largely determine the patterns of species diversity and composition [4].In contrast, stochastic theories emphasize the importance of chance colonization, random extinction and ecological drift [5].A plethora of studies supported that deterministic and stochastic processes simultaneously control the assembly of microbial communities, and their relative importance is mediated by environmental factors [6,7].However, despite knowing the structure and functions of microbial community are critical to ecosystem functioning, quantifying ecological processes controlling community composition is extremely challenging.With the rapid development of high-throughput sequencing technologies, a null modeling-based approach has been developed from the statistical perspective based on large experimental data [8].This statistical approach represents a significant advance in microbial ecology that provides insights into the assembly mechanisms in a wide variety of ecosystems [9,10].However, to our knowledge, no studies have been conducted using this approach to examine the assembly mechanisms of bacterial communities in silage ecosystem.
In this study, we evaluated the bacterial communities in oat silages successively harvested starting from the grain-filling stage through to full ripening stage.The oat was selectively harvested from southern China considering that silages produced in this area are more challenged by humid and rainy climate than temperate regions [11].To infer the assembly rules of bacterial community during silage fermentation, the statistical framework based on the null model was implemented.Furthermore, the assembly processes of abundant and rare subcommunities were, respectively, quantified, as they generally exhibit different functional traits in various ecosystems [12].This work will provide new insights into the mechanisms underlying the construction of silage bacterial community and increase our ability to harness beneficial microbiome for silage production.

Experimental material and sampling procedure
The oat was grown at Baima Teaching and Research base of Nanjing Agricultural University (31° 40′ 52″ N, 119° 5′ 28″ E), Lishui District, Nanjing, China.This area is a typical agricultural area in southern China, with an average annual temperature of 15.5-20.0°C and the average annual precipitation of 500-1000 mm.The oat was planted on 12, October 2021 in 54 experimental plots (4 m by 3 m each).All plots had the same tillage, irrigation, and fertilization practices.After 24 weeks' growth, the oat was successively harvested in 18 days starting from the grain-filling stage through to full ripening stage.The rainy days were skipped during the harvest period.On each harvest day, oat of three randomly selected plots (replicates) were harvested above 5 cm soil level and then chopped into length of 1 to 2 cm.After thorough mixing, about 500g chopped forage was sampled, and the remaining forage (approximately 500 g) was packed into polyethylene plastic bags (30 × 40 cm) and stored at room temperatures (20 to 25 °C) for 60 days.After ensiling, corresponding silage was sampled.In total, 54 raw material samples and 54 silage samples were collected for further analysis.

Experimental analysis
About 200 g of raw material and silage sample were oven dried for 48 h at 60 °C for dry matter (DM) measurements.The dried samples were then ground with a laboratory pulverizer (FW100; Taisite Instrument Co., Ltd., Tianjin, China) to pass through a 1-mm screen for total nitrogen (TN), water-soluble carbohydrate (WSC), neutral detergent fiber (NDF) and acid detergent fiber (ADF) analyses according to Dong et al. [13].The crude protein (CP) content was calculated by multiplying TN by 6.25.The DM contents were corrected with the losses of volatiles during oven drying using the equations of Gallo et al. [14].
The lactic acid bacteria (LAB), Enterobacteriaceae, aerobic bacteria, and yeasts counts were determined in raw materials according to Du et al. [15].To determine fermentation characteristics, 35 g of silage sample was blended with 60 mL distilled water and macerated for 24 h at 4 °C.The extract was filtered for pH, organic acids, and ammonia nitrogen (NH 3 -N) determinations.The pH was measured with a HANNA HI 2221 pH meter (Hanna Instruments Italia Srl, Villafranca Padovana, Italy).The NH 3 -N was determined using the phenol-hypochlorite reaction method [11].The organic acids (including lactic, acetic, propionic, and butyric acids) and ethanol were quantified using an Agilent 1260 HPLC system equipped with a refractive index detector (Carbomix ® H-NP5 column, 2.5 mM H 2 SO 4 , 0.5 mL/min).Silage quality was assessed by the index of Flieg's score with the equation: 220 + (2 × %DM-15)-40 × pH [16].According to the index, silage was considered to be very inferior when it had score of < 20; to be inferior with a score between 21 and 40; to be medium with a score between 41 and 60; to be good with a score between 61 and 80 and to be very good when it had score between 81 and 100.

DNA extraction and 16S rRNA gene sequencing analysis
DNA was extracted from the raw material and silage samples using the FastDNA SPIN Kit and the Fast-Prep Instrument (MP Biomedicals, Santa Ana, CA).The V3-V4 regions of the 16S rRNA gene of bacteria were amplified by universal primers 338F and 806R.The PCR products were purified using the AxyPrep DNA gel extraction kit (Axygen Biosciences, Union City, CA, USA).The DNA were paired-end sequenced (2 × 300 bp) on an Illumina MiSeq PE300 platform (Illumina, Inc., San Diego, CA) at Majorbio Bio-Pharm Technology Co., Ltd., Shanghai, China.Raw sequences were processed using FLASH (version 1.2.11).The QIIME (version 1.9.1) quality control process was used to eliminate sequences with quality scores below 20.Sequences with a minimum length of 200 bp were clustered into operational taxonomic units (OTUs) at a 97% similarity threshold.The assessment of alpha-diversity, including Shannon and Chao1 indexes, was carried out using QIIME (version 1.9.1).For beta-diversity analysis, principal coordinate analysis (PCoA) was performed to visualize variations in bacterial communities among samples using the UniFrac weighted-distance metric.The community structures of bacteria were analyzed at the phylum and genus levels using the Silva database (version 138) with a confidence threshold of 70%.

Assembly processes analysis
We defined OTUs as "abundant" when they had a relative abundance > 0.1%, whereas OTUs with relative abundance < 0.01% were defined as "rare" [10].All silage samples were divided into three groups (G1, G2, and G3) according to the harvest day (G1: day1-day6; G2: day7-day12; G3: day13-day18).The null model was employed to investigate the significance of deterministic and stochastic processes in generating the bacterial community structure [17].The nearest taxon index (NTI) was used to indicate whether taxa coexisting within a community are more closely related or more dispersed than would be expected by the null model.To infer the ecological processes in particular communities, we calculated the βNTI to quantify the deviation between the distribution of the βMNTD values of the observation and the βMNTD values of the null model.|βNTI|< 2 indicates that community assembly is dominated by stochastic processes.|βNTI|> 2 indicates that deterministic processes play dominant roles in community assembly.The βNTI (2) with RCbray (0.95) were combined to further determine the mechanisms of community assembly processes, such as heterogeneous selection, homogenous selection, dispersal limitation, homogenous dispersal, and undominated processes [18].

Statistical analysis
Linear regressions of chemical composition, microbial counts, or alpha-diversity estimators on harvest day were conducted and visualized by GraphPad Prism 9 (Graph-Pad Software, Inc., La Jolla, CA, USA).Statistical significance was declared at the P = 0.05 level of confidence.Principal component analysis (PCA) was performed to visualize the relationships between fermentation parameters and silages.To reveal the biomarker bacteria, the relative abundance of epiphytic and silage bacteria at genus level against harvest day were regressed using Random Forests machine learning algorithm according to Zhang et al. [19].

Raw material characteristics and fermentation quality
As shown in Fig. 1, linear increases (P < 0.001) in DM, NDF and ADF contents, aerobic bacteria, yeasts and enterobacterial counts, and linear decreases (P < 0.001) in WSC and CP contents were observed with the delay of harvest.The LAB populations fluctuated strongly during the harvest period (P = 0.1039).After 60 days of ensiling, most silages displayed "inferior" or "very inferior" fermentation qualities (Fig. S1A).There was great variation in fermentation qualities even for silages harvested on the same day.The fermentation qualities showed an overall decreasing trend with the delay of harvest (P = 0.096) (Fig. S1B).PCA showed that "very inferior" silages contained more acetic acid, propionic acid, and NH 3 -N, whereas "inferior" silages contained more butyric acid, and "very good" and "good" silages had greater lactic acid contents (Fig. S1C).

Diversities and compositions of epiphytic and silage bacterial communities
After quality control, 5,012,301 high-quality sequences were identified in the raw and silage samples, which were grouped into 1648 OTUs based on 97% sequence similarity threshold.Bacterial richness and diversity were, respectively, evaluated by Chao1 and Shannon index.The bacterial richness of epiphytic bacterial community was unaffected by the harvest day (Fig. 2A).However, epiphytic bacterial diversity linearly decreased (R 2 = 0.163, P = 0.0025) with the delay of harvest (Fig. 2B).For silage microbiota, bacterial richness linearly decreased (R 2 = 0.113, P = 0.013) (Fig. 2A), while bacterial diversity linearly increased (R 2 = 166, P = 0.002) with the delay of harvest (Fig. 2B).The PCoA and Anosim results showed that ensiling markedly altered the bacterial community (R = 0.338, P = 0.001) (Fig. 2C).Also, significant Fig. 1 Raw material characteristics across different harvest days.A-I The relationship between raw material characteristics and harvest day as revealed by simple linear analysis.Solid lines represent the linear regression models.DM: dry matter; CP: crude protein; WSC: water-soluble carbohydrates; NDF: neutral detergent fiber; ADF: acid detergent fiber; LAB: lactic acid bacteria; CFU: colony forming unit differences were observed both in epiphytic and silage bacterial communities among harvest days (epiphytic: R = 0.588, P = 0.001; silage: R = 0.423, P = 0.001) (Fig. 2D  and E).
All samples were screened and classified at phylum and genus levels.As shown in Fig. 3, the dominant epiphytic bacterial phyla were Proteobacteria and Actinobacteria.After ensiling process Firmicutes replaced Actinobacteria as the dominant phyla.Delaying harvest increased the relative abundance of Proteobacteria both in epiphytic and silage microbiota.At genus level, the dominant epiphytic genera were Rhodococcus, unclassified_f_Enterobacteriaceae, and Ochrobactrum.After ensiling the dominant bacterial genera were Lactobacillus, Pediococcus, unclassified_f_Enterobacteriaceae and Hafnia-Obesumbacterium.To identify biomarker taxa, the relative abundance of bacteria at genus level against harvest day were regressed using Random Forests machine learning algorithm (Fig. 4A and C).Totally, 13 and 6 genera were identified as biomarker taxa, respectively, in epiphytic and silage microbiota.For epiphytic microbiota (Fig. 4B), Pseudomonas, Achromobacter, Delftia, Ochrobactrum and Rhodococcus showed decreasing relative abundances with the delay of harvest, whereas the opposite trends were observed in the relative abundances of unclassified_f_Enterobacteriaceae, Erwinia, Curtobacterium.For silage microbiota (Fig. 4D), the relative abundances of Rhodococcus, Pediococcus, Enterococcus, and Bifidobacterium decreased, while the relative abundance of Hafnia-Obesumbacterium increased with the delay of harvest.The relative abundance of Lactobacillus in silages first increased and then decreased with the delay of harvest.
To reveal the relationship between measured variables and silage bacterial communities, RDA analysis and Spearman's correlation heatmap were conducted (Fig. S2).RDA analysis showed that RDA1 and RDA2 explained 35.2% of the variance in silage bacterial communities (Fig. S2A).The pH, lactic acid, acetic acid, propionic acid, DM, WSC, NH 3 -N, CP, ADF and NDF were significant predictors (P < 0.05) for the silage bacterial communities.There existed positive correlations among CP, WSC, LA and DM.These variables were negatively correlated with pH, ADF, NDF, NH 3 -N, acetic acid, and propionic acid.The Spearman's correlation heatmap showed that CP, WSC, lactic acid and DM had positive correlations with Lactobacillus and Pediococcus (P < 0.05), whereas pH, ADF, NDF, NH 3 -N, and acetic acid had positive correlations with Hafnia-Obesumbacterium and Bifidobacterium (P < 0.05) (Fig. S2B).Bacterial assembly processes and driving factors βNTI was calculated based on null model to infer the ecological processes of abundant and rare bacteria during silage fermentation (Fig. 5).The βNTI values for abundant subcommunity of G1 silages distributed between -2 and 2, suggesting that the bacterial assembly was mainly controlled by stochastic processes (Fig. 5A and C).However, with the delay of harvest, the βNTI values of abundant subcommunity increased to > 2 in G2 and G3 silages, suggesting an increased importance of deterministic processes in the community assembly.For rare subcommunity, βNTI values were mainly distributed between -2 and 2, suggesting that stochastic processes were the primary contributor to the community assembly (Fig. 5B).By combining βNTI and RCbray metrics (Fig. 5C and  D), we observed that the turnover of abundant subcommunity was mainly shaped by heterogeneous section and dispersal limitation, and rare bacteria subcommunity was mainly governed by undominated and dispersal limitation processes.In addition, the effect of heterogeneous selection on abundant subcommunity was greater in G2 and G3 silages compared with those in G1 silages (Fig. 5C).
To explore the drivers of assembly processes for bacterial taxa, linear regression analysis was applied to show the relationship between the changes of measured variables and the βNTI values of abundant bacteria (Fig. 5E-G).The βNTI values had significant correlations with the changes in DM, pH and WSC during silage fermentation (P < 0.05), suggesting that these variables are important factors in the assembly of abundant community.Among these variables, the changes in DM were positively correlated with the βNTI values (R 2 = 0.220, P < 0.001), while the changes in WSC (R 2 = 0.139, P = 0.006) and pH (R 2 = 0.099, P = 0.021) were negatively correlated with the βNTI values.Moreover, the linear regression analysis was also applied between the βNTI values and Flieg's score (Fig. 5H).We found that the βNTI values of abundant bacteria were negatively correlated with the Flieg's scores of silages (R 2 = 0.091, P = 0.026).

Raw material characteristics and fermentation qualities
Plant maturity at harvest is a crucial factor affecting the nutrient characteristics of forage [20].As expected, the DM, NDF and ADF contents increased and CP contents declined with the delay of harvest.This was mainly associated with the accumulation of structural carbohydrates and its dilution effect on CP content [21].Similarly, Zhao et al. [22] found that the increased proportions of NDF and ADF led to a relative decrease in CP content of sweet sorghum.The WSC contents also decreased with the delay of harvest.It was in consistent with Stirling et al. [23], who reported that the WSC contents in oat decreased from heading stage (49.5 g/kg DM) to water ripe stage (33.5 g/kg DM).
The population size of epiphytic bacteria on plants is largely limited by nutrients leached onto the leaf surface [24].It has been reported that the abundance of the leached nutrients from leaves tend to increase as plant matures [25].This may account for the greater aerobic bacteria, yeasts and enterobacteria populations with the delay of harvest.However, strong fluctuations of epiphytic LAB population were observed during the harvest period.It suggested that LAB population may be driven more by other temporally varying environmental factors, such as temperature, solar radiation, and ambient humidity [13,26].
Southern China is special area for silage making as hot and humid climate make the quality of silage in this area unstable and uncontrollable compared to those in northern China or temperate regions [26].In practice, silage producers have to ensile forages at relatively high moisture levels [11].To modeling the farmers' systems of practices in southern China, whole-crop oat in this study was freshly ensiled at low DM conditions (< 22% FW).The DM and WSC contents of the forage have major effects on ensiling process.According to Kung [27], undesirable bacteria, such as clostridia and enterobacteria, are easy to thrive in extremely wet silages.Our results supported this as only a few silages were identified as "very good" or "good" quality silages.Moreover, possibly due to a decrease in water-soluble carbohydrates (WSC) content, the fermentation qualities of silages showed a declining trend with the delay of harvest.This is consistent with Fang et al. [11], who reported that the change in initial WSC level could alter the fermentation pattern and affect the outcome of silage fermentation.

Diversities and compositions of epiphytic and silage bacterial communities
Epiphytic bacteria play an important role in plant growth and resistance to pathogen infection [28].Epiphytic bacterial richness displayed no statistical difference during the harvest period, suggesting that epiphytic microbiota is highly conserved during plant growth.Bacterial diversity decreased with the delay of harvest, presumably associating with the decreasing leaf ratio, since leaf contributes most to the diversity of epiphytic microbiota [29].Epiphytic microbiota structures underwent dynamic change during the harvest period.This could be due to the altered plant status including hormonal and physiological changes as matured [30].Three epiphytic genera (Enterobacter, Erwinia, Curtobacterium) exhibited increasing relative abundances with the delay of harvest.Most members of these bacteria are pathogenic to plants [31].The increases in their proportions probably reflected a decrease in plant resistance.
Generally, diverse bacterial communities are formed in field and LAB development during silage fermentation will simplify bacterial community resulting in a decline in alpha-diversity [13].Greater silage bacterial diversities with the delay of harvest suggested a decreased effect of LAB on dominating silage microbiota.The bacterial community shifted from Proteobacteria to Firmicutes after fermentation.This was closely associated with the combined stress of low pH and anaerobiosis during fermentation [32].Lactobacillus, Pediococcus, unclassified_f_Enterobacteriaceae, and Hafnia-Obesumbacterium constituted the predominated genera in the silages.Lactobacillus and Pediococcus are desirable bacteria in silage fermentation.Their flourishment has been shown to promote the establishment of acidic environment and the suppression of undesirable bacteria [33].In contrast, unclassified_f_Enterobacteriaceae and Hafnia-Obesumbacterium are the undesirable bacteria competing with LAB for limited WSC contents [11].Silages dominated by these bacteria often exhibit high pH and extensive protein degradation [11,34].This explained their positive relationships with pH and NH 3 -N contents (Fig. S2B).Delaying harvest increased the proportions of Hafnia-Obesumbacterium, explaining why the silage fermentation quality decreased with the delay of harvest.

Bacterial assembly processes and driving factors
In microbial ecology, it is widely accepted that microbiota assembly patterns in different habitats can be explained by deterministic and stochastic processes, based on niche and neutral theories, respectively [35].Null model analysis provides a way to explore whether communities are randomly assembled or non-randomly aggregated or segregated, and to identify the underlying mechanisms for microbial assembly [10,17,36].In this study, null model analysis was for the first time applied to silage ecosystem to reveal the ecological processes controlling bacterial assembly.The results showed that the assembly of rare subcommunity was primarily controlled by stochastic processes.This is likely due to the small population sizes of the rare species, which make them easily to be impacted by demographic stochasticity [37].In contrast, abundant species often occupy core niche positions and therefore they are strongly impacted by deterministic filtering [38].Our results showed that stochastic processes also played important roles in the succession of abundant subcommunity in the silages harvested in southern China.The stochastic processes consider that random changes shape microbial communities and that their fluctuations are random, including unpredictable interference, random birth and death, and dispersal probability [39].This was corroborated by the low explanatory power of measured variables in RDA analysis (Fig. S2A).It is worth noting that, deterministic processes, more specifically, heterogeneous selection, were more prominent in abundant bacteria with the delay of harvest.Generally, heterogeneous selection is determined by dynamic selection under biotic or abiotic conditions and can lead to large changes in microbial community [40].This implies that the characteristic changes during forage maturity influenced the rules governing the assembly of bacteria during silage fermentation.
To evaluate the potential drivers of producing trends in phylogenetic assembly of abundant bacteria in the oat silages, the βNTI values were correlated with the changes in measured variables during silage fermentation.The results revealed the importance of DM and WSC contents as critical factors that impact the balance between stochastic and deterministic processes in the assembly of abundant bacteria.The DM content associates with the moisture in silage, which can structure microbial communities through many indirect pathways.For example, changes in solute diffusion and water potential due to varying moisture levels contribute to distinct variations in microbial community [41].The bacterial competition for resources is regulated by the dissolved nutrient level [42].Therefore, it is not surprising that WSC content can shape the turnover of bacterial community during fermentation.In addition to DM and WSC contents, the change in pH value was also an important driver of the bacterial assembly.This could be attributed to the neutral nature of endocellular pH of most microorganisms.However, unlike our expectation, the βNTI values were negatively correlated with pH values, indicating that greater pH changes promoted stochastic assembly during silage fermentation.The disagreement between our expectation and the observation was probably explained by the fact that, besides direct effects, pH may indirectly affect the bacterial community by altering the solubility of elements (e.g., phosphorus, aluminum, and iron) [42].Bacterial species may respond differently to the direct and indirect effects of silage pH decline especially under extremely wet conditions.Similarly, studies on freshwater lakes and agricultural soils also reported the increased importance of stochasticity in acidic environments [6,43].Our results further suggested the negative relationship between Flieg's scores and βNTI values.It revealed that deterministic filter could increase the heterogeneity of community through the selection of undesirable species with stronger competitive abilities in silages.

Conclusions
This study quantified the bacterial assembly processes in oat silages harvested in southern China using null models.Significant differences in raw material characteristics were observed among harvest days.The fermentation qualities of silages tended to decrease with the delay of harvest.During silage fermentation, stochastic processes were the primary contributor to the assembly of rare subcommunity, while abundant subcommunity was controlled both by stochastic and deterministic processes.Delaying harvest increased the dominance of deterministic assembly of abundant subcommunity.The changes of three variables (DM, WSC and pH) have significant relationships with the assembly of abundant bacteria in oat silages harvested in southern China.Furthermore, significant negative correlation was found between Flieg's scores and the βNTI values.This study revealed the ecological processes controlling the bacterial assembly during silage fermentation, which provides new insights into the mechanisms underlying the construction of silage bacterial community.

Fig. 2
Fig. 2 Alpha-and beta-diversity of bacterial communities in raw materials and silages across different harvest days.A, B Linear regressions of Chao1 and Shannon indexes in raw materials and silages versus harvest day.Solid lines represent the linear regression models.C Principal component coordinate analysis (PCoA) of bacterial communities in raw material and silage samples.D PCoA of bacterial communities in raw material samples grouped by harvest day; E PCoA of bacterial communities in silage samples grouped by harvest day.Bacterial communities are compared using the unweighted UniFrac distance metric

Fig. 3 Fig. 4
Fig. 3 The bacterial community compositions in raw materials and silages across different harvest days.Bacterial community composition on phylum (A) and genus (C) level in raw materials.Bacterial community composition on phylum (B) and genus (D) level in silages

Fig. 5
Fig. 5 Bacterial assembly processes and driving factors.A The beta nearest taxon index (βNTI) values of abundant subcommunity during silage fermentation.B The beta nearest taxon index (βNTI) values of rare subcommunity during silage fermentation.C The relative contribution of deterministic and stochastic processes in the assembly of abundant subcommunity.D The relative contribution of deterministic and stochastic processes in the assembly of rare subcommunity.G1: silages harvested from day1 to day6; G2: silages harvested from day7 to day12; G3: silages harvested from day13 to day18.E-G Relationships between the βNTI values of abundant bacteria and changes in DM, WSC and pH during silage fermentation.H Relationships between the βNTI values of abundant bacteria and Flieg's scores.Solid lines in E-H represent the linear regression models.FW: fresh weight; DM: dry matter; WSC: water-soluble carbohydrates