Association of Ruminal Volatile Fatty Acids with Average Daily Gain in Growing Sheep and Goats: A Systematic Review and Meta-Regression Analysis

1. Introduction

Rumen volatile fatty acids (VFA) supply most of the energy that ruminants use, covering roughly 70%–80% of metabolizable energy (Bergman, 1990). In practice, the fermentation measures reported most often are acetic acid, propionic acid, and the A:P ratio (Russell and Hespell, 1981). These acids do not form a single uniform energy pool; instead, the VFA profile emerges from the combined action of the microbial community, the diet, and how animals are managed. Fibrous feeds typically steer fermentation toward acetate, whereas readily fermentable starch favors propionate-producing taxa such as Prevotella, Succinivibrionaceae, and Ruminococcus spp., raising the share of propionic acid (Janssen, 2010; Lin et al., 2023). Propionic acid in turn is a major substrate for hepatic gluconeogenesis: it feeds the tricarboxylic acid cycle through the propionate pathway and can contribute around 60%–74% of the carbon skeleton of glucose (Aschenbach et al., 2010). Since the propionic acid fraction reflects both fermentation type and energy partitioning, it may carry information about animal growth that other measures do not.

Average daily gain (ADG) is widely used as a growth indicator in sheep and goats (National Research Council, 2007). Calcium propionate supplementation improved lamb growth in one meta-analysis (Orzuna-Orzuna and Lara-Bueno, 2023), and goats with high daily gain showed higher ruminal propionate concentrations (Chen et al., 2024), as did high-yielding Hu sheep (Jia et al., 2026). A meta-analysis of goats likewise reported that additives such as yeast, which raise ruminal propionate and total VFA, improved ADG (Ogbuewu and Mbajiorgu, 2023); in lambs, post-weaning ADG has been connected to rumen microbiota composition and its fermentation products (Yin et al., 2023). Multi-omics work further suggests that differences in carbohydrate-degrading capacity between ruminal and hindgut microorganisms can directly affect ADG (Krone et al., 2024; Dai et al., 2021).

The literature was searched systematically in line with the PRISMA guidelines (Page et al., 2021), identifying 85 studies with 290 treatment groups. Complete data—all covariates plus the standard error of ADG—were available for 83 independent studies covering 283 treatment groups. Two analytical routes were then followed: inverse-variance-weighted multilevel meta-regression (Tempelman, 2025) as the primary analysis, and an unweighted mixed linear model as a confirmatory check. Both were used to examine how acetic acid, propionic acid, A:P, and total VFA relate to ADG in growing sheep and goats, to characterize their functional forms, and to adjust for feed intake, initial body weight, and diet composition. Incremental predictive value and the limits of applicability were assessed with the incremental heterogeneity-explanation proportion, species interaction tests, and leave-one-out cross-validation. Our hypothesis was that propionic acid adds predictive value for ADG beyond standard animal and dietary information, largely because of stable between-study differences in fermentation profiles that trace back to the biological mechanisms of the fermentation–glycogenolysis axis.

2. Materials and Methods

2.1. Data Sources and Retrieval Strategy

The systematic literature search was conducted in accordance with the PRISMA guidelines (Page et al., 2021) (search databases, search terms, and search date are detailed in Supplementary Material S1). No upper date limit was imposed; the search cutoff date was May 20, 2026. A total of 1,846 articles were identified, after deduplication, 1,275 articles remained; these underwent a preliminary screening based on titles and abstracts, eliminating 1,085 articles. Among the remaining 190 articles, full-text evaluations were performed; subsequently, 105 articles were excluded due to reasons such as failure to report ADG or VFA metrics, lack of standard deviation or standard error, experimental subjects not being sheep or goats, or duplicate publications. Ultimately, 85 studies comprising 290 treatment groups were included (48 studies involving 164 treatment groups for sheep, and 37 studies involving 126 treatment groups for goats; the literature screening process is detailed in Supplementary Material S1). A complete data sample was defined as records containing all predictive factors, covariates (DMI, IBW, CP, NDF), and the standard error (SE) for ADG; this resulted in 83 independent studies and 283 treatment groups, and Table 1A lists these 83 studies (46 sheep studies with 157 treatment groups and 37 goat studies with 126 treatment groups). For the total VFA analysis, since some studies did not report values, the analysis was based on 82 studies and 279 treatment groups. The means, standard deviations, and ranges for each predictive factor and covariate are presented in Table 1B of the article. The covariates included dry matter intake (DMI, kg/d), initial body weight (IBW, kg), dietary crude protein (CP,% DM), and neutral detergent fiber (NDF,% DM); these were independently extracted from the original manuscripts of each included study by two researchers and cross-verified; missing records were excluded according to the complete data principle. The selection of the above four factors was based on their decisive role in determining growth performance (energy intake, body size, and fiber regulation), supplemented by model comparisons incorporating stepwise addition of covariates (using AIC and heterogeneity explanation proportions) and a variance inflation factor (VIF) diagnostic (with all VIFs < 2.3) to validate their necessity and ensure non-collinearity (see Appendix A, covariate diagnostics).

2.2. Analytical Framework

The primary analysis rested on two independent frameworks. In the first, inverse variance weighting was used: each treatment group was weighted by the reciprocal of the squared standard error (SE) of its ADG mean, and multilevel meta-regression was then applied to the group means (Tempelman, 2025). Predictors were centered by their inverse-variance-weighted means, and squared terms were added to limit collinearity between X and X² and to keep the main effects easy to interpret. The random-effects structure operated on two levels, study (Study) and the within-study treatment effect, the latter being the observation-level random term that captures residual between-group heterogeneity. Treating studies as clusters, we obtained cluster-robust standard errors (CR1); P-values under this framework are therefore cluster-robust, while within-model P-values are reported for reference. The second framework was an unweighted linear mixed model (LMM) in which Study entered as a random intercept and effect size (ID) was nested within Study as a random effect; marginal R² was derived with the Nakagawa–Schielzeth method. The inverse-variance framework was designated the primary (formal) analysis and supplied the principal results, and the unweighted LMM served as a confirmatory check that the findings did not hinge on the weighting choice. Because the LMM's random structure already accounts for within-study correlation, within-model P-values are reported for this framework. Heterogeneity explained is summarized by three complementary proportions (see Figure 1): the proportion for the model with predictors only (M1; the “(X only)” value in Figure 1) relative to the empty model; the incremental proportion of the full model relative to the covariates-only model (incremental R², reported untruncated so that negative values are kept); and the proportion of the full model relative to the empty model. The output also includes τ² for each model (spanning both study and treatment components), the residual heterogeneity Q-test, and the overall I². These three measures index different benchmarks—completeness, incrementality, and overall heterogeneity—so they should not be compared directly.

$Y_{\mathit{ij}}=A_0+A_i+A_{1i}{\mathit{Goat}}_{\mathit{ij}}+A_{2i}{\mathit{Pro}}_{\mathit{ij}}+A_{3i}{\mathit{DMI}}_{\mathit{ij}}+A_{4i}{\mathit{IBW}}_{\mathit{ij}}+A_{5i}{\mathit{CP}}_{\mathit{ij}}+A_{6i}{\mathit{NDF}}_{\mathit{ij}}+m_{\mathit{ij}}$

Where: Yij is the average daily gain (ADG, kg/d) of the jth treatment group in the ith study; $A_0$ is the overall intercept; $A_i$ is the random intercept for study i (random effect); $A_{1i}$ is the species-effect coefficient, with "Goat" an indicator variable taking 1 for goats and 0 for sheep; $A_{2i}$, $A_{3i}$, $A_{4i}$, $A_{5i}$, and $A_{6i}$ are the fixed-effect coefficients for propionic acid concentration (Pro, mmol/L), dry matter intake (DMI, kg/d), initial body weight (IBW, kg), dietary crude protein (CP, % DM), and neutral detergent fiber (NDF, % DM), respectively; and mᵢⱼ is the residual term.

$\mathit{ADG}(\mathrm{kg}/d)=0.0349-0.0055·\mathit{Goat}+0.00219·\mathit{Pro}+0.1320·\mathit{DMI}-0.00131·\mathit{IBW}+0.00136·\mathit{CP}-0.00115·\mathrm{NDF}$

Standard errors are shown in parentheses: intercept 0.0349 (0.0356); Goat = −0.0055 (0.0133); Pro = 0.00219 (0.00037); DMI = 0.1320 (0.0144); IBW = −0.00131 (0.00090); CP = 0.00136 (0.00085); NDF = −0.00115 (0.00048). n = 273 treatment groups and 80 studies; the marginal R² (R²_m) = 0.658.

2.3. Automatic Selection of Function Forms

To determine whether a linear or quadratic functional form should be adopted to assess the association between each predictor (acetic acid, propionic acid, A:P, total VFA) and ADG, after controlling for covariates, linear and quadratic models were fitted separately for each indicator (with the quadratic term included as a centered variable to reduce multicollinearity); multiple pieces of evidence were then synthesized to decide on the appropriate functional form: statistical significance of the quadratic coefficient (P < 0.01); the incremental heterogeneity explanation proportion of the quadratic model relative to the linear model under either the inverse-variance framework (Figure 1) or the mixed linear framework (Figure 2), or the relative increase in the marginal R² (≥10%); and the degree of improvement in the Akaike Information Criterion (AIC) (a decrease ≥4). The quadratic form was selected when it outperformed the linear form for an indicator; otherwise, the linear form was retained to avoid over-parameterization and extrapolation beyond observed data ranges. The results showed that, under the inverse-variance framework, the functional forms for A:P and total VFA were identified as quadratic, while those for acetic acid and propionic acid were linear; under the mixed linear framework, all indicators exhibited linear relationships.

2.4. Species Interaction Test

To explore the sources of heterogeneity further, the primary analysis was extended by fitting an expanded multilevel meta-regression model with average daily gain (ADG) as the outcome and species, dry matter intake (DMI), initial body weight (IBW), dietary crude protein (CP), and neutral detergent fiber (NDF) as covariates. To separate within-study from between-study sources of treatment effects, the treatment-level covariate (X) was decomposed into a within-study component (X_within—the deviation of each treatment group's X from the mean X of its own study) and a between-study component (X_between—the mean X across all studies); main effects and their species interactions were estimated separately for the two components. Weighting (inverse variance) and the random structure (studies as random intercepts, effect ID nested within each study) matched the primary analysis. For the interaction term, slope differences between sheep and goats were tested at both the within-study and between-study levels, with multiple interaction tests corrected using the BH method. When an interaction was not significant, a contrast approach estimated within-study and between-study slopes for each species with their 95% confidence intervals (CIs), and the overlap ratio among the ranges of original X, within-study X, and between-study X quantified the degree of data overlap between the two species.

2.5. Leave-One-Sample Cross-Validation (LOSO)

For each predictor, a baseline model (DMI + IBW + CP + NDF; species also included in the full sample) and a baseline-plus-indicator model were estimated. Leave-one-out cross-validation was then applied: each iteration omitted one study, refit the model to the remaining studies, and predicted all treatment groups of the omitted study. Evaluation proceeded along three lines: overall population prediction (sheep and goats), a sheep-specific equation, and a goat-specific equation. The cross-validation R² (CV-R²), the change in CV-R² after adding the indicator (ΔCV-R²), the RMSE, and its percentage reduction were computed. Paired bootstrap resampling (1000 repetitions) with the study as the sampling unit supplied 95% confidence intervals (CI) for the reductions in ΔCV-R² and RMSE, taken as the 2.5% and 97.5% percentiles. A 95% CI for ΔCV-R² that excludes zero together with a lower RMSE denotes clear incremental predictive value; if only the point estimate improves while the 95% CI crosses zero, predictive performance improves but the increment remains uncertain.

2.6. Sensitivity Analysis

The SE Monte Carlo perturbation approach drew 500 independent random values for the standard errors of the predictors, DMI, and IBW under a normality assumption, re-estimated the interaction model each time, and reported whether the perturbed percentile intervals excluded zero and whether the signs matched the primary findings. The full dataset and the complete-data subset were compared by repeatedly fitting the primary model to the entire dataset with the data-based approach, checking whether the direction, significance, and incremental R² stayed qualitatively unchanged. Robustness was further assessed by comparing coefficients, P-values, and signs across four model variants: inverse variance weighting, a model with only study-level random effects, precision truncation (95th percentile), and the unweighted mixed-effects model. Publication bias was appraised with a funnel plot combined with Egger's linear regression test, in which standardized effect sizes were regressed on the ADG standard errors of the studies to test whether the intercept deviated from zero. Egger's test P-values for all four indicators were < 0.001, hinting at possible small-sample or publication bias; association magnitudes may accordingly be inflated, and the results should be interpreted with this caveat (see Supplementary Material S1). Individual study risk of bias was not formally assessed with a dedicated tool (e.g., ROBINS-I); instead, reporting bias was evaluated via funnel plots and Egger’s test, and the robustness of inferences was examined through multiple sensitivity analyses (Section 2.6). The review protocol was not registered in PROSPERO.

2.7. Statistical Software, Significance Definition, and Data Availability

All analyses were performed in R (version 4.5.1). Multilevel meta-regression models were fitted with the rma.mv function of the metafor package (version 4.6.0); linear mixed models (LMMs) were fitted with the nlme package (version 3.1-166) and the lme4 package (version 1.1-35.5); data management and variable transformations relied mainly on the dplyr package (version 1.1.4). The significance level was set at α = 0.05; the association between the four VFA indicators and ADG was treated as a pre-specified primary hypothesis; P-values were used to assess the robustness of the main effects; and multiple comparisons for species interactions were corrected with the BH method. Analysis scripts, diagnostic plots, and all intermediate results are provided in Appendix A and Supplementary Material S1.

3. Results

3.1. Sample Characteristics

The complete data sample comprised 83 studies and 283 treatment groups (46 studies with 157 treatment groups in sheep and 37 studies with 126 treatment groups in goats); the total VFA dataset included 82 studies and 279 treatment groups. The baseline characteristics and literature sources of the included studies are presented in Table 1A. Within the complete data sample, the mean ADG was 0.170 ± 0.081 (SD) kg/d, the mean acetic acid was 48.30 ± 22.70 mmol/L, the mean propionic acid was 17.95 ± 9.83 mmol/L, the mean A:P was 3.06 ± 1.15, and the mean total VFA was 77.27 ± 36.77 mmol/L; the mean ADG values were 0.213 kg/d for sheep and 0.117 kg/d for goats. The means, standard deviations, and ranges of all variables are provided in Table 1B. A parallel analysis using feed conversion rate (FCR, i.e., feed intake-to-gain ratio) as a secondary endpoint is available in Supplementary Materials S2–S5.

3.2. The Correlation Between Acetic Acid and ADG

Acetic acid showed a significant positive association with ADG (Table 2A: β = 0.00088; 95% CI: 0.00031–0.00144; cluster-robust P = 0.003; total model heterogeneity explained = 0.583), and this relationship was also significant in the mixed linear framework (Table 2B: β = 0.00070; 95% CI: 0.00034–0.00106; P < 0.001; marginal R² = 0.610). The incremental heterogeneity explained by acetic acid in Table 2A was 0.017, the smallest among the four indicators. At the species level, sheep showed significant associations in both frameworks (Table 2A: β = 0.00098; robust P = 0.013; Table 2B: β = 0.00066; P = 0.004), whereas goats did not in either (Table 2A: β = 0.00052; robust P = 0.205; Table 2B: β = 0.00063; P = 0.065). See Figure 1 and Figure 2.

3.3. Association Between Propionic Acid and ADG

Propionic acid showed a significant positive correlation with ADG (Table 2A: β = 0.00188, 95% CI: 0.00095–0.00280; cluster-robust P < 0.001; overall model heterogeneity explanation proportion = 0.621); the mixed linear framework was also significant (Table 2B: β = 0.00213, 95% CI: 0.00139–0.00287, P < 0.001; marginal R² = 0.637). In Table 2A, the incremental heterogeneity explanation proportion for propionic acid was 0.107, the highest among the four indicators. Within the species subgroup, both frameworks for goat propionic acid yielded significant results (Table 2B: β = 0.00279, P < 0.001; Table 2A: β = 0.00227, robust P < 0.001); for sheep propionic acid, the mixed linear framework yielded a significant result (Table 2B: β = 0.00157, P = 0.002), while the inverse variance framework also yielded a significant result (Table 2A: β = 0.00148, robust P = 0.035).

3.4. The Relationship Between A:P and ADG

The evidence strength for A:P is inconsistent across the two frameworks: in the inverse-variance framework, it does not achieve cluster-robust significance (Table 2A: β = –0.00820; 95% CI: –0.01720 to 0.00080; robust P = 0.074), whereas in the mixed linear framework, it exhibits weak negative significance (Table 2B: β = –0.00769; 95% CI: –0.01384 to –0.00153; P = 0.015). The incremental heterogeneity explanation proportion is 0.071 (Table 2A). The contrast—where concentration-based indicators show a positive association while proportion-based indicators show a negative association—is consistent with the exploratory findings from the between-study decomposition analysis; however, the inconsistency in significance across the two frameworks suggests that the A:P result is not robust and should therefore not be used as an independent predictive indicator.

3.5. Correlation Between Total VFA and ADG

The total VFA and ADG exhibited a significant positive correlation (Table 2A: β = 0.00088, 95% CI: 0.00042–0.00135; cluster-robust P < 0.001; the incremental heterogeneity explanation proportion for total VFA reached 0.083). The mixed linear framework also showed a significant result (Table 2B: β = 0.00127, 95% CI: 0.00080–0.00174, P < 0.001; marginal R² = 0.642). At the species level, both frameworks yielded significant results for sheep (Table 2A: β = 0.00108, robust P = 0.002; Table 2B: β = 0.00157, P < 0.001) and for goats (Table 2A: β = 0.00055, robust P = 0.033; Table 2B: β = 0.00148, P = 0.039). For details, see Figure 1 and Figure 2.

3.6. Species Interaction Test-

In the extended model, neither the Species × X (Within) nor the Species × X (Between) interaction term reached statistical significance (P > 0.05). Species-specific slopes showed that for sheep acetic acid both the intra-study and inter-study slopes were significant (P = 0.005 and P = 0.032), whereas for goat acetic acid neither slope was significant. For goat propionic acid both the intra-study and inter-study slopes were significant (P = 0.001 and P = 0.031), while for sheep propionic acid only the inter-study slope was significant (P = 0.008). These findings describe a point-estimate pattern in which acetic acid is the more sensitive indicator for sheep and propionic acid for goats; nevertheless, the differences were not significant enough to confirm an interaction effect.

3.7. Incremental Predictive Value of Leave-One-Study-Out Cross-Validation (Table 3)

When the overall equation was used for prediction, adding propionic acid lifted CV-R² from 0.393 to 0.518 (ΔCV-R² = 0.124; 95% CI: 0.012–0.326; excludes zero) and cut RMSE by 10.84%, so both metrics point to clear incremental predictive value. Total VFA (ΔCV-R² = 0.070) and A:P (ΔCV-R² = 0.073) improved the point estimates, but because their 95% CIs crossed zero the improvements were not statistically significant (Figure 3). For sheep, propionic acid raised CV-R² to 0.485 (ΔCV-R² = 0.056; 95% CI crosses zero) and lowered RMSE by 5.03%, an improvement whose incremental benefit is still uncertain. For goats, the same addition moved CV-R² from −0.599 to −0.157 (ΔCV-R² = 0.441; 95% CI: 0.035–1.337; excludes zero) and reduced RMSE by 14.92%, which we read as clear incremental predictive value even though the goat baseline was inherently weak. Within both species equations, point estimates for each indicator were best when propionic acid was included (ΔCV-R² of 0.046 and 0.204, respectively), yet both 95% CIs crossed zero; these results should be treated as exploratory. Figure 4 plots observed against predicted values after propionic acid was added to the overall equation.

Table 3. Incremental predictive value of leave-one-study-out (LOSO) cross-validation.

Validation Scenario	Number of groups	Baseline CV-R²	propionic acid-model CV-R²	ΔCV-R² (propionic acid vs. baseline, 95% CI)	RMSE reduction (%)	propionic acid P
Overall→Overall	279	0.393	0.518	0.124(0.012,0.326)	10.84	< 0.001
Overall→Sheep	153	0.429	0.485	0.056(−0.056,0.165)	5.03	< 0.001
Overall→Goat	126	−0.599	−0.157	0.441(0.035,1.337)	14.92	< 0.001
Sheep→Sheep	153	0.384	0.431	0.046(−0.052,0.160)	3.82	< 0.001
Goat→Goat	126	−0.250	−0.047	0.204(−0.061,0.547)	8.51	< 0.001

In the column headings, A→B denotes the model trained on dataset A and validated on dataset B; the baseline model contained DMI, IBW, CP and NDF (plus species for the total sample); LOSO provides internal validation.;ΔCV-R² = CV-R² (with propionic acid) − baseline CV-R² (without propionic acid); the 95% CI was derived from study-level paired bootstrap (1,000 resamples; 2.5% and 97.5% percentiles); bold indicates that the 95% CI excludes 0.;N denotes the number of treatment groups predicted in each validation direction; P for propionic acid is the P-value of the propionic acid coefficient in the LMM that included propionic acid.;A negative baseline CV-R² indicates that the baseline model predicted worse than simply using the overall mean, reflecting the weak baseline predictive ability for goats.

Table 3. Incremental predictive value of leave-one-study-out (LOSO) cross-validation.

Validation Scenario	Number of groups	Baseline CV-R²	propionic acid-model CV-R²	ΔCV-R² (propionic acid vs. baseline, 95% CI)	RMSE reduction (%)	propionic acid P
Overall→Overall	279	0.393	0.518	0.124(0.012,0.326)	10.84	< 0.001
Overall→Sheep	153	0.429	0.485	0.056(−0.056,0.165)	5.03	< 0.001
Overall→Goat	126	−0.599	−0.157	0.441(0.035,1.337)	14.92	< 0.001
Sheep→Sheep	153	0.384	0.431	0.046(−0.052,0.160)	3.82	< 0.001
Goat→Goat	126	−0.250	−0.047	0.204(−0.061,0.547)	8.51	< 0.001

In the column headings, A→B denotes the model trained on dataset A and validated on dataset B; the baseline model contained DMI, IBW, CP and NDF (plus species for the total sample); LOSO provides internal validation.;ΔCV-R² = CV-R² (with propionic acid) − baseline CV-R² (without propionic acid); the 95% CI was derived from study-level paired bootstrap (1,000 resamples; 2.5% and 97.5% percentiles); bold indicates that the 95% CI excludes 0.;N denotes the number of treatment groups predicted in each validation direction; P for propionic acid is the P-value of the propionic acid coefficient in the LMM that included propionic acid.;A negative baseline CV-R² indicates that the baseline model predicted worse than simply using the overall mean, reflecting the weak baseline predictive ability for goats.

3.8. Sensitivity Analysis

Comparing the complete-data and full-data results, the direction of effects was identical across all four indicators, and the qualitative significance conclusions were also consistent. The Monte Carlo perturbation test (500 replicates) and the multi-inference robustness check support the main conclusion; the propionic acid interaction and the predictive increment remained directionally stable under SE perturbations. When only the two inference variants—the study-level random-effects model and the 95th-percentile precision ceiling—were considered, the directions and significance of acetic acid, propionic acid, and total VFA were largely consistent; however, neither acetic acid nor propionic acid was statistically significant (see Supplementary Material S1, Sensitivity Table).

4. Discussion

We combined 83 studies and 283 treatment groups under two complementary frameworks—inverse-variance-weighted meta-regression and multilevel mixed-effects models—to ask how total VFA, acetic acid, propionic acid, and the A:P ratio relate to average daily gain (ADG) in growing sheep and goats, and whether the VFA profile adds predictive value on top of conventional covariates. The findings separate association from prediction rather cleanly. On the association side, acetic acid, propionic acid, and total VFA were positively related to ADG in both frameworks, whereas A:P was inconsistent. On the prediction side, only propionic acid produced a clear incremental cross-study gain (ΔCV-R² = 0.124; RMSE reduced by 10.84%), and the species-specific point estimates never translated into significant interaction effects. The rest of the discussion develops these results at three levels—mechanism, species, and application.

4.1. Principal Analysis Results and Mechanisms

4.1.1. Separation of Correlation Strength and Predictive Increment

Across both statistical frameworks, total VFA, acetic acid, and propionic acid were significantly and positively correlated with ADG, and after adjustment for DMI, IBW, CP, NDF, and species the direction of the effects stayed consistent with similar magnitudes (Table 2A and Table 2B). This cross-method agreement is a first sign of reproducibility. For the three indicators, the full model explained 0.625, 0.583, and 0.621 of the between-study heterogeneity, respectively. Turning to incrementality, the VFA profile added only modest heterogeneity beyond conventional covariates, though all values were positive: 0.017 for acetic acid, 0.071 for A:P, 0.107 for propionic acid, and 0.083 for total VFA, with propionic acid again the largest. LOSO reinforced this picture, showing the clearest cross-study predictive gain for propionic acid (ΔCV-R² = 0.124; 95% CI: 0.012–0.326; RMSE reduced by 10.84%). In short, propionic acid outperformed the other variables on robust association, incremental R², and external prediction. The concentration-based indicators (acetic acid, total VFA), by contrast, tracked the scale–mass dimension that covaries with feed intake, body size, and diet composition, so their incremental explanatory power was limited. Predictive value should therefore be judged by incremental R² or external validation rather than by single-variable correlations, consistent with earlier studies that found VFA relate to growth indicators yet add little predictive gain. For A:P, evidence strength differed between the two frameworks and the point estimates were small, making the result unstable. The functional form also diverged between frameworks (quadratic in Table 2A; linear in Table 2B), so a linear model is the most straightforward and appropriate choice in both cases. Its interpretation and application should be confined to the range of the observed data.

4.1.2. Mechanistic Reasoning Behind the Relationship Between Fermentation Products and Tissue Growth

For the statistical structure described above to carry meaning, every link in the chain must rest on established biochemical and physiological evidence rather than on mere correlation with phenotypes. The framework therefore runs through four interconnected processes: microbial production, epithelial absorption, portal vein and hepatic metabolism, and endocrine and tissue synthesis. Production comes first: fermentation kinetics set the metabolic profile, with each sugar molecule degraded into acetic, propionic, and butyric acids in characteristic proportions, and the propionic acid pathway competing with methanogenesis for hydrogen (Janssen, 2010; Morgavi et al., 2010). Starch-based diets favor taxa such as Succinivibrionaceae and Prevotella, which raises propionic acid levels (Lin et al., 2023; Krone et al., 2024); elevated propionic acid is at once an ecological outcome of intense fermentation and a key precursor for hepatic gluconeogenesis. Absorption comes next. VFA cross the epithelium by passive diffusion or through transporters such as MCT1/MCT4, NHE, and NBCe1; most butyric acid is oxidized in the epithelium to β-hydroxybutyrate, whereas propionic and acetic acids reach the portal vein in relatively higher proportions (Aschenbach et al., 2011; Rathert-Williams et al., 2023). Epithelial morphology and transporter expression adapt to fermentation load as well (Penner et al., 2011; Na and Guan, 2022). One point deserves emphasis: rumen fluid concentration is not the same as substrate flux reaching the liver. Absolute concentration depends heavily on rumen fluid volume, saliva secretion, and sampling time, while the molar ratio is largely insensitive to dilution—an anatomical reason why propionic acid prediction is more robust and why concentration-based indicators are more exposed to measurement noise (Fortina et al., 2022). Metabolically, propionic acid is converted by propionyl-CoA carboxylase (PCC) to methylmalonyl-CoA and then isomerized to succinyl-CoA by methylmalonyl-CoA mutase (MCM) with vitamin B12 as coenzyme; succinyl-CoA is later turned into oxaloacetate through the reverse tricarboxylic acid cycle (rather than by direct catalysis via PCK1/PEPCK), moved into the cytosol by the mitochondrial malate/aspartate shuttle, and finally converted by phosphoenolpyruvate carboxykinase (PCK1) to phosphoenolpyruvate (PEP), which enters gluconeogenesis to form glucose and supplies 60%–74% of the carbon skeleton (Aschenbach et al., 2010). The same pathway upregulates hepatic gluconeogenic enzymes and has insulinotropic effects (Wang et al., 2023; Pang et al., 2023). Chen et al. (2026) showed in a lamb gavage experiment that higher ruminal propionic acid upregulates hepatic PC and G6PC and boosts hepatic glucose output (Chen et al., 2026); the glucose–insulin–IGF-1 axis then drives nitrogen deposition and protein synthesis. Because higher propionic acid means more glucose precursors, propionic acid sits at a metabolic hub in the positive causal chain. Feedback closes the loop: propionic acid oxidation also acts as a satiety signal (hepatic oxidation theory), so the ADG effect of propionic acid operates within a homeostatic loop of energy supply and appetite regulation, with the marginal effect flattening as concentration rises rather than growing without bound (Allen, 2014); very high fermentation loads may even trigger subacute ruminal acidosis that lowers feed intake and growth (Golder and Lean, 2024; Elmhadi et al., 2022). Table 2A shows no significant quadratic term for propionic acid, so within the observed range the relationship looks predominantly linear. Still, the mechanisms above hint that endogenous effects of very high propionic acid may decelerate. This logic explains both the strong concentration-based correlations we observed and the prediction of stable structural configurations, the stability being owed to the measurement configuration rather than to the strength of the biological associations.

4.2. Species Differences and Cross-Study Predictive Value

The species-specific slopes follow a point-estimation pattern in which acetic acid is the more sensitive indicator for sheep and propionic acid for goats, yet none of the Species × X interaction terms reached statistical significance. We therefore treat these results as exploratory rather than fitting separate models per species. Variation in microbial community structure among sheep breeds may also weaken the stability of the species-specific signal (Xiang et al., 2022). The pattern fits with the host's dual control over fermentation dynamics: goats graze more often, have a relatively larger rumen volume, and show more flexible retention times, so the effective fermentation intensity and profile differ between the species even on identical diets (Krone et al., 2024). The host genome is a major determinant of the core microbiota and markedly shapes the rumen microbial community (Lin et al., 2023; Krone et al., 2024; Dai et al., 2021), and the microbial response to dietary change is fairly host-specific (Henderson et al., 2015; Weimer, 2015). The non-significant interaction is consistent with Salah et al. (2015), who reported that intake responses in sheep, goats, and cattle run parallel but differ in intercepts, systematic between-species differences appearing as baseline rather than slope differences (Salah et al., 2015). A unified framework with species-level effects is therefore preferable to separate models per species. On a point-estimation basis, LOSO indicated that propionic acid was the optimal additive component in every predictive scenario: the gain in overall prediction was clear (ΔCV-R² = 0.124; 95% CI excludes zero; RMSE reduced by 10.84%), and the improvement was largest for goats (ΔCV-R² = 0.441; RMSE reduced by 14.92%) despite a weak baseline; in all other scenarios the 95% CI crossed zero. It is worth noting that the incremental gain was greatest for targets with weaker baselines: when conventional covariates fail to capture the underlying pattern, configuration indicators may carry more marginal value, and such an inverse relationship between predicted gain and baseline performance is visible only through externalized leave-one-out cross-validation.

4.3. Comparison with Previous Studies and Robustness

A number of independent findings support the positive correlation we observed between propionic acid and ADG. Lambs given calcium propionate grew better in a meta-analysis (Orzuna-Orzuna and Lara-Bueno, 2023); goats with high daily weight gain carried higher ruminal propionic acid levels (Chen et al., 2024); high-yielding Hu sheep had elevated ruminal propionate (Jia et al., 2026); yeast-based supplements that raised ruminal propionate and total VFA improved ADG in goats (Ogbuewu and Mbajiorgu, 2023); post-weaning ADG in lambs was tied to ruminal microbiota composition (Yin et al., 2023); and an oral administration experiment confirmed that increased ruminal propionate upregulates key hepatic gluconeogenic enzymes (Chen et al., 2026). The reverse pattern was also observed: Velázquez-Cruz et al. (2024) found that low doses (10 g/kg DM) of calcium or sodium propionate failed to improve growth in fattening lambs (Velázquez-Cruz et al., 2024), matching our view that the effects depend on dose and physiological stage—the growth-promoting action of propionate is not unconditional, and the ADG response to propionic acid should be read as the net result of the fermentation–metabolism chain rather than of the additive alone. Earlier studies generally reported effects through correlation or regression P-values, while few provided incremental R², robust inference, and leave-one-out validation together. The present study addresses that gap by reporting the limited explanation of incremental heterogeneity transparently and by using LOSO to quantify the cross-study predictive benefit of propionic acid directly. Robustness checks support the same conclusions: the inverse-variance framework (Table 2A) and the mixed linear model framework (Table 2B) agreed on the significant positive correlation of acetic acid, propionic acid, and total VFA, and when the A:P direction was treated as non-robust, neither the complete-versus-full dataset comparison, nor the SE Monte Carlo simulation (500 replicates), nor the multi-simulation inference changed the direction or qualitative conclusion (Result 3.8), indicating that the primary findings do not hinge on specific data-cleaning or statistical decisions.

4.4. Application Significance and Limitations

The concentration- and composition-based approaches described above—designed to capture fermentation intensity and compositional stability in prediction—point to a differentiated application strategy. Total VFA and acetic acid work as monitoring indicators of fermentation intensity, suited to tracking fermentation scale in ration formulation and livestock management. Propionic acid emerges as the more robust compositional indicator across studies (overall equation prediction ΔCV-R² = 0.124, RMSE reduced by 10.84%), which agrees with evidence that propionic acid production improves feed efficiency and growth performance (Jia et al., 2026; Ogbuewu and Mbajiorgu, 2023; Yin et al., 2023; Chen et al., 2026). A:P, in turn, should not be used alone to predict growth performance. The Total–propionic acid model can serve as an empirical estimation model built on treatment-group means and checked through leave-one-out internal validation; it estimates mean performance across treatment groups under similar experimental conditions, but it does not support precise individual prediction and calls for independent external validation. Some design limitations should be kept in mind when interpreting these results. Because the findings rest on treatment-group means, publication bias is a real risk; funnel plots and Egger's regression flagged publication bias for all four indicators (Egger's P-values < 0.001; see Supplementary Material S1), which may have inflated the association strength. Heterogeneity arising from breed, feeding practices, sampling times, and VFA determination methods was not fully accounted for (Fortina et al., 2022), and the total VFA dataset covered only 279 treatment groups across 82 studies, with the functional-form choice sensitive to the analytic framework. Leave-one-out validation is internal validation, so the model's predictive power for baseline values in goats is relatively weak and does not justify developing or relying on goat-specific equations. Taken together, these limitations argue that the quantitative inferences should be read at the level of average effects across studies, with caution when extrapolating to individual animals.

5. Conclusion

This study used a multilevel modeling and cross-validation framework to assess how well the VFA profile predicts ADG. Total VFA, acetic acid, and propionic acid were robustly associated with ADG in both analytical frameworks; beyond conventional covariate adjustment, however, the VFA profile added only limited incremental explanatory power for inter-study heterogeneity. Among the components, propionic acid showed the highest incremental explanatory value and the most stable additional contribution at the cross-study prediction level, whereas A:P was neither robust nor independently informative.

Point estimates differed between sheep and goats, but the differences were not statistically significant at the interaction level. We therefore recommend incorporating propionic acid into a growth prediction model based on treatment-group means, with priority given to validation using independent external data; application of such models should be restricted to estimates derived from treatment-group means under conditionally similar settings. Future research should combine mechanistic validation with a richer set of host covariates to deepen our understanding of VFA-mediated growth regulation.