Estimating energy concentrations in wooded pastures of NW Spain using empirical models that relate observed metabolizable energy to measured nutritional attributes

sustainability Article Estimating Energy Concentrations in Wooded Pastures of NW Spain Using Empirical Models That Relate Observed Metabolizable Energy to Measured Nutritional Attributes María Pilar González-Hernández 1,* and Juan Gabriel Álvarez-González 2   Citation: González-Hernández, M.P.; Álvarez-González, J.G. Estimating Energy Concentrations in Wooded Pastures of NW Spain Using Empirical Models That Relate Observed Metabolizable Energy to Measured Nutritional Attributes. Sustainability 2021,13, 13581. https://doi.org/10.3390/su132413581 Academic Editors: Paris Fokaides and Chia-Lin Chang Received: 23 September 2021 Accepted: 6 December 2021 Published: 8 December 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Department of Crop Production and Projects of Engineering, Higher Polytechnic Engineering School, University of Santiago de Compostela, 27002 Lugo, Spain 2Unit of Sustainable Environmental and Forest Management, Department of Agroforestry Engineering, Higher Polytechnic Engineering School, University of Santiago de Compostela, 27002 Lugo, Spain; [email protected] *Correspondence: pilar[email protected] Abstract: Wooded pastures serve as a traditional source of forage in Europe, where forest grazing is valued as an efficient tool for maintaining the diversity of semi-natural habitats. In a forest grazing setting with diverse diet composition, assessing the energy content of animal diets can be a difficult task because of its dependency on digestibility measures. In the present study, prediction equations of metabolizable energy (ME) were obtained performing stepwise regression with data (n= 297; 44 plant species) on nutritional attributes (Acid Detergent Fiber, lignin, silica, dry matter, crude protein, in vitro organic matter digestibility) from 20 representative stands of Atlantic dry heathlands and pedunculate oak woodlands. The results showed that the prediction accuracy of ME is reduced when the general model (R 2 = 0.64) is applied, as opposed to the use of the specific prediction equations for each vegetation type (R 2 = 0.61, 0.66, 0.71 for oak woodlands; R 2 = 0.70 heather-gorse dominated heathlands, R 2 = 0.41 continental heathlands). The general model tends to overestimate the ME concentrations in heaths with respect to the observed ME values obtained from IVOMD as a sole predictor, and this divergence could be corrected by applying the specific prediction equations obtained for each vegetation type. Although the use of prediction equations by season would improve accuracy in the case of a Winter scenario, using the general model as opposed to the prediction equations for Spring, Summer or Fall would represent a much smaller loss of accuracy. Keywords: acid detergent fiber; Atlantic dry heathlands; forest grazing; in vitro organic matter digestibility; pedunculate oak woodlands; silvopasture 1. Introduction Grazing is one of the influential and sustainable management interventions for shaping and maintaining semi-natural habitats and promoting greater biodiversity and multifunctionality in forest understory and heathlands [ 1 – 4 ]. Horse grazing, for example, can decrease gorse dominance and promote heath composition in Atlantic heathlands [ 4 – 6 ], which are prized for their biodiversity, aesthetic and cultural value, and supported by the European Union (EU) within the High Nature Value farming framework. Grazing oak woodlands is included as a type of silvopasture that has been practiced throughout the ages across Europe and has gained growing recognition in the last decades [ 7 – 12 ]. Several studies have also pointed out the importance of livestock grazing as a management strategy that can improve rural sustainability while effectively controlling the accumulation of flammable woody vegetation derived from rural depopulation and the abandonment of management practices in the heathlands of NW Spain [ 5 , 6 , 13 , 14 ], thus contributing to the mitigation of greenhouse gas emissions by reducing the risk of fire [ 15 ]. Furthermore, in projected climate change scenarios, a stable forage supply in hot, dry years would be Sustainability 2021,13, 13581. https://doi.org/10.3390/su132413581 https://www.mdpi.com/journal/sustainability Sustainability 2021,13, 13581 2 of 22 provided only by extensive and moderate farming, which allows the development of an insulating tree cover within grazed pastures [16]. Thus, various reasons motivate the growing interest in shrubs as fodder in spite of their limited digestibility, particularly in heathlands [ 17 ]. Traditional rustic cattle breeds in the north-west of the Iberian Peninsula (e.g., Vianesa, Frieiresa, Cachena, Maronesa or Arouquesa) and goats are adapted to browse on plants considered to be of “low” nutritional value that are present in the understory of forest stands and in heathlands [ 12 , 18 , 19 ]. Such browsing can have multiple benefits in animal nutrition, including increased feed intake, digestibility and rumen fermentation, and reduced methane emission from ruminants [ 20 ]. Extensive grazing with these rustic and traditional livestock breeds is recognized for its important environmental, social, cultural, market and public value, and many of these breeds have the protected geographical status from the European Commission because they have suffered a drastic population decline in past decades [ 21 ]. The current social demand for organic meat obtained from free-range livestock management that involves browse with low inputs of supplemental feed resources adds up to this general growing interest in wooded pastures. On the other hand, decisions regarding optimum grazing densities are often connected with the nutritional composition of pastures, as well as their intake and renewal rates [ 22 – 24 ], and those factors can be crucial in avoiding poor achievement of management targets that can arise from inappropriate grazing. The nutritional profile of food selected by herbivores is often positively related with moisture content, crude protein (CP) and metabolizable energy (ME), which are among the most important constituents in predicting food selection [ 25 ]. There is a high correlation between ME and fiber content due to the fact that it decreases the digestibility on which ME depends. A low time retention of high fiber diets in the animal digestive tract increases the dry matter voluntary intake and translates into a more rapid pasture consumption, affecting the food-resource renewal rates as well as the dynamic of changes in vegetation composition [ 22 , 24 ]. For all these reasons, data on nutrient availability is required to minimize the environmental impact of animal production and range management. Unfortunately, accurate measurement of the ME values of any given diet requires the execution of in vivo digestibility trials with animals, and sometimes the measurement of the Gross Energy (GE) of feed, feces, urine and methane. Most such measurements can only be made in research centers specially prepared to do so, as well as being costly; consequently, to solve these difficulties, some approaches to estimated energy concentrations from nutrient composition or in vitro digestibility measurements have been proposed [ 26 , 27 ]. Because of the relative constancy of both the GE of digestible organic matter for most forages other than silage and the proportions of energy lost in urine and methane, highly significant correlations are normally found between ME and organic matter digestibility (OMD), such that approximate conversions of in vitro OMD (IVOMD) values to ME values may be made [ 28 ]. Nonetheless, the determination of IVOMD has greater limitations than other chemical parameters such as acid detergent fiber (ADF) or its components (lignin, cellulose and silica), and crude protein (CP); both of these are related to the rate of digestibility and are more readily available to the livestock manager. In extensive forest grazing scenarios, where the diet is composed of a great variety of plant species characterized by different nutritional attributes, the estimation of ME concentrations becomes a complicated task; in addition, few studies have reported on the value of browse digestibility for heathlands and pedunculate oak woodlands [ 17 , 29 – 32 ]. All these challenges have motivated a growing interest in the prediction of the energy of diets by establishing relationships between observed ME and nutritional parameters that are easily analyzed [ 33 – 35 ]. Conversions of IVMOD data into ME values [ 28 ] provide an estimated ME that is possible to relate to measured nutritional parameters using empirical models, thus facilitating a rapid prediction of ME concentrations adapted to a particular vegetation type. Such predictive equations are very useful when a marginal error in the accuracy of estimates is tolerable, although it means examining a large number of samples widely distributed in time and space. Sustainability 2021,13, 13581 3 of 22 The aim of the current study was to obtain prediction equations of ME using a wide range of data on nutritional attributes (IVOMD, ADF, cellulose, lignin, silica, dry matter and crude protein) from various stands that represent prototypes of European Atlantic heathlands and pedunculate oak woodlands. We hypothesized that if the existing correlations between the measured nutritional characteristics and observed ME could result in obtaining significant empirical linear models for such vegetation types, then such prediction equations could be broadly applicable in stands with similar vegetation depending on their available data on chemical composition, thus facilitating grazing management decisions that are based on the ME requirements of the animals. 2. Materials and Methods 2.1. Study Area Twenty forest stands located in Galicia, north-west Spain (42 ◦ 45 0 N, 7 ◦ 51 0 W), a region of the Iberian Peninsula biogeographically assigned to the European Atlantic province [ 36 ], were selected for the study (Figure 1). These represent three archetypal botanical associations of pedunculate oak (Quercus robur) woodlands(10sites)and two typesofheathlands: Atlantic heathgorse heathlands (8) and continental dry (2) heathlands. The climate of the area is Atlantic, with mild, wet winters at the coast, and colder inland. Mean annual temperatures in oak woodlands and Atlantic heathlands (AH) were between 9.7–14.2 ◦ C, decreasing to 7.3–8.3 ◦ C in continental heathlands (CH). Annual precipitation in the sites ranged from 1419–2100 mm. Elevation of the study sites varied from 60–800 m for oak woodlands and Atlantic heathlands, whereas continental heathland sites reached 1120–1280 m. Soils are mainly acidic, having developed from granite, slate or schists. The prevalent shrubs that comprise oak woodland understory vary according to the following three botanical associations [ 37 ]: Blechno spicant-Quercetum roboris (O- Bl), Myrtillo-Quercetum roboris (O-My) and Rusco aculeati-Quercetum roboris (O-Ru). The dominant shrubs differed in the two types of heathlands, with the presence of Atlantic components such as Ulex spp., Daboecia cantabrica, Calluna vulgaris and various Erica spp. in the heath-gorse type distinctive from a less diverse shrub layer in the continental type. More details about the key species considered for the present study within each vegetation type can be seen in Tables A1 and A2 (Appendix A). Figure 1. Location of the 20 forest stands selected for the study. O-Bl : (1) Boimorto, (2) Cerceda, (3) Cerqueiras, (4) Fragavella, (5) Monfero. O-My : (6) A Lama, (7) Lovios, (8) Xestoso. O-Ru : (9) A Rua, (10) Lourizán. AH : (11) Arnuide, (12) Cerqueiras, (13) Cerponzóns, (14) Coto do Muiño, (15) Fragavella, (16) Oia, (17) Pouzos, (18) Silleda. CH: (19) Eirexa, (20) Rubias. Sustainability 2021,13, 13581 4 of 22 2.2. Nutritional Characterization of Vegetation For the nutritional analysis, 44 plant species were selected according to their abundance in each forest site and their palatability for ruminants (e.g., roe deer, goat, beef cattle, sheep and horse). Their occurrence by vegetation type was as follows: 20 such species were present in O-Bl, 17 in O-My, 12 in O-Ru, 16 in Atlantic heathlands and 5 in continental heathlands. Plant samples for each species were composed of apical portions no more than 15 cm in length and 0.5 cm in diameter, and were clipped randomly from selected individuals before being combined into a single sample for analysis. Samples were collected quarterly during 2 years for most of the species. For nutritional characterization, acid detergent fiber (ADF), permanganate lignin, cellulose and silica contents were obtained with sequential detergent analysis [ 38 ]. Crude protein was determined by analyzing the nitrogen content using the micro-Kjeldahl method, with the results being multiplied by 6.25. In vitro organic matter digestibility (IVOMD) was assessed using the technique of Tilley and Terry (1963) as modified by Alexander (1969) [ 39 , 40 ]. More details about the methods of sampling and analysis can be found in González-Hernández and Silva-Pando (1999) [ 17 ]. Additionally, the IVOMD values were used to estimate observed ME using the equation reported by Barber et al. (1984) and recommended by AFRC (1993) [ 28 ]: ME (MJ/kg DM) = 0.0157 × IVOMD (g/kg DM). ME concentrations obtained from applying this equation will be referred to in the text as observed ME concentrations or ME equations based on IVOMD as a sole predictor. With the purpose of being a reference for the nature and range of the data from which the ME equations were obtained, the nutritional characteristics by plant species are summarized in Appendix A. We used the coefficient of variation (CV) to refer to the relationship between the size of the mean and the variability of each nutritional parameter. Values of observed ME for the same plant species were grouped into those cases where such species occurred in more than one type of oak woodlands, or was in the two types of heathlands, and when certain homogeneity in their ME values was observed (i.e., low CV). In the field of agricultural sciences, CV is broadly accepted when it is ≤ 20%. A scale of reference is often as follows: ≤ 10 is very good, 10–20 is good, 20–30 is acceptable and CV ≥ 30 is not acceptable. Because of the limited records of some forbs and grasses in our data set, we decided to pool them together. 2.3. Data Analysis to Obtain Prediction Equations of ME To analyze the possible linear relationships between ME and nutritional characteristics, Pearson correlation coefficients were calculated and the trends of pairs of ME data with each possible dependent variable were visualized. In view of these trends, it was decided to also include the logarithmic transformation of nutritional characteristics as possible dependent variables. Cellulose was excluded as an independent variable because it is a parameter rarely consider in the nutritional reports from laboratory analysis. It is also included as a fraction of the parameter ADF, and its occurrence as a good predictor in attempts of ME prediction equations is not common. Empirical linear models that relate the observed ME to the measured nutritional characteristics were fitted for each one of the 5 vegetation types for each season, and for the pooled data using the SAS/STAT ® software [ 41 ]. Construction of the models involved the following two steps: (i) selection of the best set of independent variables and (ii) evaluation and cross-validation of the fitted models by graphical and numerical analysis of the residuals. To select the best set of independent variables, the stepwise variable selection method [ 42 ], combined with the values of the Pearson correlation coefficients and an understanding of the fitting process, was used. Multicollinearity is a problem that could be associated with the proposed models that violate the fundamental least squares assumption of independence and equal distribution of errors with zero mean and constant variance. Multicollinearity refers to the existence of high intercorrelations among the independent variables in multiple linear or nonlinear regression analysis, because some of the variables represent or measure similar Sustainability 2021,13, 13581 5 of 22 phenomena. Although the least squares estimates of regression coefficients remain unbiased and consistent under the presence of multicollinearity, this problem may seriously affect the standard errors of the coefficients, invalidating confidence intervals [ 43 ]. Thus, models with low multicollinearity should be selected whenever possible [ 44 ]. To evaluate the presence of multicollinearity among variables in the models analyzed, the variance inflation factors (VIF) of all the independent variables were calculated, and values up to 5 were accepted [45]: VIF =1 1−R2 i where R 2i is the multiple determination coefficient obtained when the ith independent variable Xiis regressed against all of the remaining independent variables in the model. 2.4. Model Evaluation and Validation Four statistical criteria obtained from the residuals (e i ) were examined to test the accuracy and precision of ME estimates: bias ( E ), mean absolute error (| E |), root mean square error (RMSE) and the adjusted coefficient of determination (R 2adj ). These expressions can by summarized as follows: Bias =E= n ∑ i=1 (yi−ˆ yi) n MAE =E= n ∑ i=1 |yi−ˆ yi| n RMSE =s∑n i=1(yi−ˆ yi)2 n−1 R2 adj =1−∑n i=1(yi−ˆ yi)2·(n−1) ∑n i=1(yi−y)2·(n−p) where ei=yi− _ yi and y i , ˆ yi and y are the measured, predicted and average values of the dependent variables, respectively; nis the total number of observations used to fit the model and pis the number of model parameters. Because the quality of the fit does not necessarily reflect the quality of the prediction, cross-validation of each model was carried out, and the residual for one observation was estimated by fitting the model without that observation. MEF, equivalent to the R 2 of the fitting phase, was then calculated from the residuals from the cross-validation. Finally, the estimates obtained with the general model (fitted to all data) when applied to each of the five vegetation types and to the different stations are compared with the estimates obtained with the models fitted specifically for those situations in order to estimate the error of using the general model when it is not possible to discriminate in the database between vegetation types or seasons. 3. Results 3.1. Nutritional Characteristics, Digestibility and ME (Based on IVDOM as Sole Predictor) The chemical composition and digestibility of the vegetation types and plants studied here are summarized in Appendix A. A high variation in observed ME concentrations in continental and Atlantic heathlands was found (CV = 35.9 and 45.6%, respectively), suggesting a large range in the observed ME concentration for the main plants included in these vegetation types. Observed ME values for continental heathlands ranged between 1.7–7.2 MJ/kg DM, whereas ME concentrations of plants in Atlantic heathlands extended to a wider 1.8–12.2 MJ/kg DM interval. The lowest ME concentrations in heathlands were found in heaths, followed by Halimium spp. and leguminous shrubs. On the other hand, the variation in ME values in understory plants of O-Bl and O-Ru oak woodlands was smaller than in heathlands (CV 28.7, Sustainability 2021,13, 13581 6 of 22 32.8%), and ranged from 2.2–12.6 and 3.6–11.2 MJ/kg DM, respectively. Only the degree of variation in ME values in O-My (CV = 38.3%; 2.8–13.6 MJ/kg DM) exceeded the one found for continental heathlands. ME concentrations of oak woodland understory is expected to be higher than heathlands because of the superior ME contribution from plants like Hedera helix,Ilex aquifolium,Lonicera periclymenum, leguminous brooms and Frangula alnus. Forbs and grasses, when present, would contribute considerably to increasing the ME concentrations in all vegetation types (see Tables A1–A4 in Appendix A). The variation in ADF, lignin and cellulose contents for each vegetation type was smaller than the ME oscillations, especially for heathlands (CV between 9.6–17%). The CP contents and CP variation found in continental heathlands (50.9 and 86.8 g/kg DM, CV = 14%) were more minor than in oak woodlands and Atlantic heathlands (CP minimum values around 50 g/kg DM and maximum around 160–205 g/kg DM, CV between 27–30%). 3.2. Relationships between Observed Metabolizable Energy and Nutritional Parameters Pearson correlation coefficients were conducted to evaluate the relationships between nutritional composition and observed ME values (Figure 2). The correlation matrices for each vegetation type revealed that there was a negative and highly significant relationship between ME and total fiber (ADF) or its fraction cellulose (r values ranging mostly from − 0.50 to − 0.67; p< 0.001). Conversely, no significant correlation was found between ME and lignin. The relationship of silica and ME was also significant, though it was much weaker in general (r = − 0.24 to − 0.41; p< 0.01; p< 0.05). Dry matter contents of plants were also highly significant and negatively correlated with ME (r values varying from − 0.48 to − 0.70; p< 0.001). CP, the only parameter positively correlated with ME, exhibited a weaker relationship in oak woodlands, especially in O-RU (r = 0.30; p< 0.05), although this relationship was highly significant in Atlantic heathlands (r = 0.77, p< 0.001). Overall, the correlation coefficients analyses performed for each vegetation type showed similar relationships between ME and nutritional parameters, except for continental heathlands. For this vegetation type, silica and dry matter contents were not significantly related to ME. Figure 2. Pearson correlation coefficients between ME (MJ/kg DM) and nutritional parameters (g/kg DM) for pooled data (n= 297). *** indicates a statistical significance < 0.001 and * indicates a statistical significance < 0.05. 3.3. Prediction Equations of ME from Nutritional Parameters Prediction equations for ME concentration in oak woodland understory and heathlands developed in the present study are shown in Table 1. The stepwise regression analysis indicated that ADF and log DM were the best single predictors of ME in oak woodland understory. The accuracy of ME prediction for this vegetation type improved by Sustainability 2021,13, 13581 7 of 22 adding lignin and silica (as well as their natural logarithm transformations) as predictor variables, and final R 2 reached values between 0.61–0.71. In O-Ru and O-My, the final R 2 value with additional predictors rather than one single predictor rose by over 0.36 and 0.27 units, respectively. Conversely, however, the magnitude of increase for O-Bl was only 0.07 units higher. On the other hand, equations from multi-step regression analysis for heathlands revealed CP as the best single predictor for ME (Table 1). The accuracy of prediction for ME could be improved by adding ADF and silica in the latter equations (R 2 increased from 0.60–0.70 and from 0.30–0.41 in Atlantic and continental heathlands, respectively). For all vegetation types, each parameter included on the finally selected models had a highly significant effect on the relationship (p< 0.05), except the intercept of the model fitted for CH (Equation (5b)); therefore, the model was refitted without intercept (Equation (5c)). Stepwise regression analysis of the data uniquely by season (not taking into account the particular vegetation type) indicated that the best single predictors of ME were DM for nutritional data obtained in Spring and Summer, and ADF for those from Fall and Winter (Table 2). The accuracy of prediction increased by adding CP and fractions of fiber (lignin and silica) in the latter equations. Those equations showed the lowest R 2 values for Summer (0.57) and better adjustments for Spring and Winter (R 2 0.71 and 0.73, respectively). Again, for all seasons, each parameter had a highly significant effect on the relationship (p< 0.05). The study of residuals allowed the comparison of the degree of change in the accuracy of ME estimations among the prediction equations obtained in this study. The use of the prediction equation obtained from total pooled data (i.e., general model or Equation (6d)) to estimate the ME in each vegetation type showed RMSE values higher than those found when operating with specific prediction equations for each vegetation type (Table 3). The highest increase of error was shown for Atlantic heathlands (29%) and the lowest for continental heathlands (2%). The rise in error when using the Equation (6d) to estimate ME for each type of oak woodland (i.e., instead of its specific equation) was around 11%. Similarly, the use of prediction equations obtained by season to estimate ME showed a lower RMSE than when the general model was applied. The expected loss of accuracy from using the general model instead of the prediction equations obtained for each season is expected to be especially high for estimates of ME concentrations in the winter (almost 25% RMSE increment). For ME estimates of vegetation in the Fall and Summer, the use of a general model versus a seasonal specific model would represent a loss of accuracy of 4%, which becomes only 1.8% in Spring. 3.4. Comparison of Predicted ME from Equations Developed in the Present Study and Using IVOMD as a Sole Predictor The range of predicted ME concentrations in oak woodlands narrowed with respect to observed ME concentrations (3.4–12 MJ/kg DM for O-Bl, 4.3–11 MJ/kg DM for O-Ru), except for O-My (2.3–13.8 MJ/kg DM). In Atlantic heathlands the minimum and maximum values increased for predicted ME (2.7–12.9 MJ/kg DM), and for continental heathlands the amplitude of range was reduced to values between 2.6 and 6.5 MJ/kg DM. The results of graphing the observed ME against predicted ME (equations developed in the present study) revealed the degree of disparity expected in the ME values obtained when using the prediction equations as opposed to one based on its IVOMD as a sole predictor (Figures 3–5). Sustainability 2021,13, 13581 8 of 22 Table 1. Stepwise prediction of ME from the chemical composition of understory in three types of oak woodlands (10 sites), in Atlantic heath-gorse (8 sites), Continental dry (2 sites) heathlands and for the pooled data (all vegetation types). All variables left in the selected model are significant at the 0.05 level. Data Set 1Equation Code Equations R2MEF RMSE (MJ/kg DM) O-Bl (n= 69) 1a ME = 37.06541 −5.21387 ln DM 0.54 1.5143 1b ME = 34.82655 −0.00418 ADF − 4.41128 ln DM 0.57 1.4637 1c ME = 26.75274 −0.00613 ADF + 1.69057 ln Lignin − 4.31669 ln DM 0.61 0.55 1.3919 O-My (n= 54) 2a ME = 13.87868 −0.01348 ADF 0.44 1.8033 2b ME = 29.98207 −0.01009 ADF − 3.07073 ln DM 0.64 1.4462 2c ME = − 14.41271 − 0.02772 ADF + 8.72691 ln ADF −3.19398 ln DM 0.66 1.4035 2d ME = − 22.59736 − 0.03427 ADF + 10.22191 ln ADF + 0.32157 ln Silica −2.93126 ln DM 0.68 1.3614 2e ME = 29.63567 −0.01175 ADF − 0.02373 Silica + 0.68538 ln Silica − 3.02345 ln DM 0.71 0.64 1.3047 O-Ru (n= 48) 3a ME = 12.86392 −0.01074 ADF 0.30 1.6879 3b ME = 36.28577 −0.01121 ADF − 4.00149 ln DM 0.59 1.2957 3c ME = 31.50200 −0.01309 ADF + 0.01434 Lignin −3.38594 ln DM 0.66 0.60 1.1748 AH (n= 93) 4a ME = −0.92486 + 0.06869 CP 0.60 1.3010 4b ME = 3.17646 −0.00534 ADF + 0.06102 CP 0.65 1.2087 4c ME = 2.82355 −0.00410 ADF − 0.01075 Silica + 0.06015 CP 0.68 1.1615 4d ME = 27.95771 −0.00281 ADF − 0.01101 Silica + 0.14590 CP − 7.55436 ln CP 0.70 0.66 1.1227 CH (n= 33) 5a ME = −19.64457 + 5.70339 ln CP 0.30 1.3237 5b ME = − 10.01671 − 0.00888 ADF + 4.84651 ln CP 0.42 1.2007 5c ME = −0.01142 ADF +2.87803 ln CP 0.41 0.36 1.2117 Pooled Data (n= 297) 6a ME = 35.17172 −4.96649 ln DM 0.44 1.7823 6b ME = 32.85982 −0.00929 ADF − 3.62794 ln DM 0.60 1.5008 6c ME = 30.69233 −0.01053 ADF + 0.00988 Lignin −3.40714 ln DM 0.63 1.4540 6d ME = 18.88889 −0.01008 ADF + 0.00951 Lignin −2.57947 ln DM + 1.51553 ln CP 0.64 0.63 1.4194 ME: metabolizable energy, ADF: acid detergent fiber, DM: dry matter, CP: crude protein. 1 Pooled data = all data for the 5 vegetation types Sustainability 2021,13, 13581 9 of 22 Table 2. Stepwise prediction of ME from chemical composition of pooled data (understory of oak woodlands (10 sites) and representative Atlantic heath-gorse (8 sites) and Continental dry (2 sites) heathlands by season. All variables left in the selected model are significant at the 0.05 level. Season Equations R2MEF RMSE (MJ/kg DM) Spring (n= 73) ME = 36.78326 −5.14401 ln DM 0.61 1.7801 ME = 36.25662 − 0.01835 Silica − 4.95451 ln DM 0.66 1.6458 ME = 35.44320 + 0.00560 Lignin − 0.01853 Silica −4.96216 ln DM 0.68 1.6188 ME = 50.51671 −3.73678 ln ADF + 0.01452 Lignin − 0.00518 Silica − 3.83576 ln DM 0.72 1.5131 ME = 53.08717 −4.35111 ln ADF + 0.01594 Lignin −3.68131 ln DM 0.71 0.68 1.5201 Summer (n= 104) ME = 33.20858 −4.67851 ln DM 0.34 1.8197 ME = 61.78940 −5.05363 ln ADF − 4.06425 ln DM 0.54 1.5094 ME = 49.65541 −4.81496 ln ADF − 3.38040 ln DM + 1.48910 ln CP 0.57 0.53 1.4729 Fall (n= 80) ME = 74.90378 −10.76743 ln ADF 0.51 1.3518 ME = 82.29102 −8.63380 ln ADF − 3.54878 ln DM 0.63 1.1571 ME = 73.04114 −6.48335 ln ADF − 0.40105 ln Silica −4.10694 ln DM 0.66 0.63 1.1152 Winter (n= 40) ME = 83.41177 −12.12545 ln ADF 0.65 1.2702 ME = 70.48888 −10.53676 ln ADF + 0.03610 CP 0.73 0.69 1.1154 ME: metabolizable energy, ADF: acid detergent fiber, DM: dry matter, CP: crude protein. Table 3. Analysis of the accuracy of prediction for ME when using the general model (derived from pooled data) as opposed to using the models obtained by vegetation type or by season. Prediction Equation Bias MAE RMSE R2RMSE% % RMSE Increment General model applied by O-Bl 0.2472 1.2880 1.5467 0.52 22.8 11.1 vegetation type O-My −0.0849 1.1139 1.4549 0.63 23.1 11.5 O-Ru 0.1633 1.0092 1.3032 0.58 18.5 10.9 AH −0.2480 1.2076 1.4510 0.50 32.0 29.2 CH 0.0858 0.9542 1.2247 0.40 28.1 1.1 Total 0.0003 1.1490 1.4194 0.64 24.6 15.3 General model applied by Spring 0.0820 1.2212 1.5474 0.70 23.0 1.8 season Summer −0.2694 1.2489 1.5314 0.53 27.1 4.0 Fall 0.1563 0.9896 1.1624 0.64 21.7 4.2 Winter 0.2405 1.0764 1.3929 0.58 27.4 24.9 Total 0.0003 1.1490 1.4194 0.64 24.6 5.4 Model by vegetation typeO-Bl 0.0022 1.1173 1.3919 0.61 20.5 O-My 0.0011 1.0137 1.3047 0.71 20.7 O-Ru 0.0029 0.9203 1.1748 0.66 16.7 AH −0.0017 0.8768 1.1227 0.70 24.8 CH −0.0012 0.9393 1.2117 0.41 27.8 Total 0.0005 0.9713 1.2315 0.73 21.4 Model by season Spring 0.0006 1.2257 1.5201 0.71 22.6 Summer 0.0000 1.2400 1.4729 0.57 26.0 Fall 0.0000 0.9236 1.1152 0.66 20.8 Winter 0.0003 0.8427 1.1154 0.73 22.0 Total 0.0002 1.0978 1.3464 0.68 23.4 Sustainability 2021,13, 13581 16 of 22 tends to overestimate the ME concentrations in heaths when compared to the ME values obtained from IVOMD as a predictor (i.e., AFRC ME estimation), and this divergence could be corrected to some extent by applying specific prediction equations for each vegetation type. Although the use of prediction equations by season would improve accuracy in the case of a Winter scenario, using the general model as opposed to the prediction equations for Spring, Summer or Fall would represent a much smaller loss of accuracy. Finally, some of the prediction equations here obtained included silica as a predictor, which has been reported as a reliable digestibility internal marker, especially when the content in plants is at least 7.0 g/kg DM. The range of silica values provided for different species in the present study will also serve as a reference to decide on the adequacy of applying equations that include silica for predicting ME in similar vegetation types. Author Contributions: Conceptualization, M.P.G.-H. and J.G.Á.-G.; methodology, M.P.G.-H. and J.G.Á.-G.; software, J.G.Á.-G.; validation, J.G.Á.-G.; formal analysis, J.G.Á.-G.; investigation, M.P.G.- H.; writing—original draft preparation, M.P.G.-H.; writing—review and editing, M.P.G.-H.; visualization, M.P.G.-H. and J.G.Á.-G.; supervision, M.P.G.-H. and J.G.Á.-G. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Data sharing is not applicable to this article. Acknowledgments: The previous laboratory and field work, from which the current study was later based on, was funded by the Rural Development Agency of Galician Government (Xunta de Galicia, Spain). Conflicts of Interest: The authors declare no conflict of interest. Sustainability 2021,13, 13581 17 of 22 Appendix A Table A1. Nutritional composition, in vitro digestibility and metabolizable energy (g/kg DM, g/kg, MJ/kg DM) of plant species in oak woodlands of NW Spain. For each species, first value corresponds to mean value, followed by standard deviation and coefficient of variation (in bold). VEGETATION TYPE SPECIES ADF LIGNIN CELLULOSE SILICA DM IVOMD CP ME Shrubs O-Bl, O-My, O-Ru Hedera helix 431.88 158.38 264.86 8.62 305.95 621.53 82.58 9.76 n= 17 71.93 40.44 38.23 4.93 38.52 34.12 15.15 0.54 16.66 25.53 14.43 57.17 12.59 5.49 18.34 5.49 O-Bl, O-My, O-Ru Rubus sp. 560.21 168.51 336.76 54.94 355.96 308.75 102.89 4.85 n= 16 93.87 38.87 59.48 39.07 62.14 55.80 21.83 0.88 16.76 23.06 17.66 71.12 17.46 18.07 21.22 18.07 O-Bl, O-My Erica arborea 607.07 160.23 433.12 13.68 443.95 238.62 78.31 3.75 n= 13 68.18 28.04 73.61 12.59 73.15 49.03 11.10 0.77 11.23 17.50 16.99 92.04 16.48 20.55 14.18 20.55 O-Bl, O-My Ilex aquifolium 570.89 203.71 356.53 10.96 377.62 477.79 83.87 7.50 n= 14 61.98 22.14 47.97 6.16 41.08 39.44 12.64 0.62 10.86 10.87 13.45 56.23 10.88 8.26 15.08 8.26 O-Bl, O-My Vaccinium myrtillus 589.88 176.12 387.45 26.28 410.92 353.80 75.86 5.55 n= 10 94.46 32.90 69.99 27.88 42.06 59.78 8.70 0.94 16.01 18.68 18.06 106.08 10.24 16.90 11.47 16.90 O-Bl, O-Ru Lonicera periclymenum 513.52 181.93 314.17 17.40 290.13 503.75 87.10 7.91 n= 12 61.18 18.12 57.94 13.16 56.18 91.97 31.67 1.44 11.91 9.96 18.44 75.61 19.36 18.26 36.36 18.26 O-My Genista florida 586.20 169.50 406.00 10.70 460.00 490.00 179.80 7.69 n= 1 O-Ru Cytisus striatus 419.9 69.6 348.2 2.1 336 500 153.1 7.85 n= 1 O-Ru Ruscus aculeatus 545.03 129.30 404.78 10.93 427.17 375.67 93.78 5.90 n= 6 79.83 43.44 44.62 6.43 74.58 28.00 20.49 0.44 14.65 33.59 11.02 58.83 17.46 7.45 21.85 7.45 Trees O-Bl, O-My, O-Ru Pyrus cordata 567.73 150.12 377.03 40.56 415.93 362.40 74.22 5.69 n= 10 122.25 36.30 93.42 32.09 68.20 71.38 13.43 1.12 21.53 24.18 24.78 79.11 16.40 19.70 18.10 19.70 O-Bl, O-My Quercus robur 640.70 218.75 403.90 18.05 460.00 320.00 90.54 5.02 n= 2 40.87 9.55 23.19 8.13 42.43 0.00 2.61 0.00 6.38 4.36 5.74 45.05 9.22 0.00 2.88 0.00 O-Bl, O-Ru Castanea sativa 532.42 162.26 316.76 53.40 378.40 423.60 105.37 6.65 n= 5 100.72 22.55 67.43 44.23 52.73 45.62 14.60 0.72 18.92 13.90 21.29 82.83 13.94 10.77 13.86 10.77 O-Bl Betula alba 657.55 152.20 400.95 104.35 370.50 315.00 108.75 4.95 n= 2 56.36 49.78 44.90 61.16 43.13 14.14 7.42 0.22 8.57 32.71 11.20 58.61 11.64 4.49 6.83 4.49 O-My Fagus sylvatica 690.54 160.22 442.90 87.33 450.56 227.00 78.38 3.56 n= 5 74.40 64.49 86.31 53.67 50.90 11.45 18.18 0.18 10.77 40.25 19.49 61.46 11.30 5.04 23.19 5.04 O-Ru Frangula alnus 513.30 187.84 309.56 15.88 325.20 563.40 110.33 8.85 n= 5 81.42 25.89 63.23 11.32 68.79 33.73 29.10 0.53 15.86 13.78 20.43 71.27 21.15 5.99 26.38 5.99 O-Ru Laurus nobilis 608.64 146.88 425.02 36.72 435.40 346.20 71.65 5.44 n= 5 46.01 22.79 57.72 20.53 65.60 45.69 15.31 0.72 7.56 15.52 13.58 55.90 15.07 13.20 21.37 13.20 Forbs O-Bl, O-My, O-Ru Teucrium scorodonia 593.90 156.36 356.82 80.70 324.36 424.15 91.44 6.66 n= 13 80.14 47.04 52.19 47.54 60.25 78.72 23.05 1.24 13.49 30.08 14.62 58.91 18.57 18.56 25.21 18.56 O-Bl, O-My Asphodelus albus 369.10 128.10 236.30 4.70 118.76 789.80 179.58 12.40 n= 5 102.80 67.49 45.06 1.89 30.99 19.52 43.78 0.31 27.85 52.68 19.07 40.17 26.09 2.47 24.38 2.47 O-Bl, O-My Other forbs* 463.48 146.00 276.35 41.04 193.18 598.00 145.55 9.39 n= 4 131.41 49.88 74.97 48.10 54.77 136.77 53.81 2.15 28.35 34.16 27.13 117.20 28.35 22.87 36.97 22.87 Grasses O-Bl, O-My, O-Ru Grasses** 655.15 170.52 419.23 65.39 324.30 423.71 100.14 6.65 n= 21 111.95 50.14 70.66 34.32 113.39 93.91 28.61 1.47 17.09 29.40 16.85 52.48 34.96 22.16 28.57 22.16 ADF: Acid Detergent Fiber, CP: crude Protein, IVOMD: In vitro Organic Matter Digestibility, ME Metabolizable Energy. Vegetation type refers to the following oak woodlands: Blechno spicanti-Quercetum roboris (O-Bl), Myrtillo- Quercetum roboris (O-My) and Rusco aculeati- Quercetum roboris (O-Ru). Same plant species occurring in more than one type of oak woodlands were grouped when their ME data showed similar (low SD and CV). In those cases, this is indicated in the vegetation type column. Other forbs* include the following species: Anemone nemorosa,Omphalodes lusitanica,Potentilla erecta,Polygonatum odoratum. Grasses** were pooled and included: Carex remota,Agrostis capillaris, A. curtisii,Brachypodium sylvaticum,Deschampsia flexuosa,Holcus mollis,Pseudarrenatherum longifolium,Luzula henriquesi. Sustainability 2021,13, 13581 18 of 22 Table A2. Nutritional composition, in vitro digestibility and metabolizable energy (g/kg DM, g/kg, MJ/kg DM) of plant species in Atlantic heath-gorse and continental dry heathlands (AH, CH) of NW Spain. For each species, first value corresponds to mean value, followed by standard deviation and coefficient of variation (in bold). VEGETATION TYPE SPECIES ADF LIGNIN CELLULOSE SILICA DM IVOMD CP ME Shrubs AH, CH Erica australis 718.55 166.05 526.98 25.50 512.28 162.08 60.11 2.54 n= 12 28.40 19.31 29.32 9.54 69.94 36.11 7.53 0.57 3.95 11.63 5.56 37.41 13.65 22.28 12.53 22.28 AH, CH Halimium lasianthum 716.90 171.38 457.31 88.19 406.07 207.25 75.58 3.25 n= 16 55.57 41.48 39.68 38.44 83.75 34.12 12.38 0.54 7.75 24.20 8.68 43.58 20.62 16.46 16.38 16.46 AH, CH Pterospartum tridentatum 590.91 183.32 397.76 9.81 541.32 349.00 74.47 5.48 n= 9 35.02 17.36 31.84 7.10 48.14 32.92 8.09 0.52 5.93 9.47 8.00 72.32 8.89 9.43 10.86 9.43 AH Calluna vulgaris 652.47 132.91 452.75 66.80 459.17 245.91 62.10 3.86 n= 11 99.68 26.21 59.36 33.95 61.63 32.85 8.32 0.52 15.28 19.72 13.11 50.83 13.42 13.36 13.40 13.36 AH Daboecia cantabrica 644.73 128.99 428.35 87.38 403.63 245.38 85.06 3.85 n= 8 113.87 42.47 67.90 41.25 56.65 53.43 15.32 0.84 17.66 32.93 15.85 47.22 14.04 21.78 18.02 21.78 AH *Erica spp. 657.54 155.24 476.52 25.72 456.55 205.21 65.42 3.22 n= 29 72.63 25.19 57.66 15.68 73.66 46.32 8.36 0.73 11.05 16.23 12.10 60.96 16.13 22.57 12.78 22.57 AH ** Ulex spp. 677.14 144.55 522.38 10.18 434.90 392.30 92.73 6.16 n= 20 61.63 31.95 46.48 7.90 77.22 55.38 15.73 0.87 9.10 22.10 8.90 77.55 17.76 14.12 16.96 14.12 CH Halimium umbellatum 630.3 194.3 414.38 21.58 448.88 304.25 68.50 4.78 n= 4 35.81 8.46 31.91 11.41 43.91 35.91 5.47 0.56 5.68 4.35 7.70 52.90 9.78 11.80 7.99 11.80 Forbs AH Asphodelus albus 417.77 143.70 267.20 6.88 172.83 767.33 151.30 12.05 n= 3 127.68 76.73 46.06 6.08 8.98 13.20 16.73 0.21 30.56 53.39 17.24 88.35 5.20 1.72 11.06 1.72 AH Potentilla erecta 542.44 186.64 322.46 33.34 319.38 415.60 101.62 6.52 n= 5 126.42 72.53 49.64 22.30 80.65 88.38 26.53 1.39 23.31 38.86 15.39 66.88 25.25 21.26 26.11 21.26 Grasses AH, CH Agrostis curtisii 700.28 188.81 452.55 58.88 407.10 392.75 78.25 6.17 n= 8 33.31 17.56 29.17 17.44 142.54 49.94 23.04 0.78 4.76 9.30 6.45 29.62 35.01 12.72 29.45 12.72 ADF: Acid Detergent Fiber, CP: crude Protein, IVOMD: In vitro Organic Matter Digestibility, ME Metabolizable Energy. Same plant species occurring in more than one type of heathland were grouped when their ME data showed similar (low SD and CV). In those cases, this is indicated in the vegetation type column. * Erica spp. (n): Erica ciliaris (7), E. cinerea (12), E. mackaiana (2), E. umbellata (8). ** Ulex spp. (n): Ulex europaeus (9), U. gallii (1), U. minor (10). Sustainability 2021,13, 13581 19 of 22 Table A3. Nutritional composition, in vitro digestibility and metabolizable energy (g/kg DM, g/kg, MJ/kg DM) for each vegetation type. VEGETATION TYPE ADF LIGNIN CELLULOSE SILICA DM IVOMD CP ME O-Bl Mean 579.27 173.01 369.54 36.79 348.16 431.70 94.24 6.78 Min 276 56.70 192.50 2.60 117.40 140.00 49.38 2.20 Max 755.3 253.00 536.30 147.60 585.00 803.00 208.50 12.61 SD 110.47 42.41 84.64 37.80 96.93 141.60 32.52 2.22 CV 19.07 24.51 22.90 102.74 27.84 32.80 34.51 32.80 O-My Mean 562.98 164.81 359.03 39.11 372.00 400.72 88.12 6.29 Min 258.80 74.40 162.50 0.20 80.00 180.00 54.94 2.83 Max 755.30 260.20 519.40 150.80 533.70 806.00 205.50 12.65 SD 118.31 39.79 85.10 41.14 99.94 153.32 30.43 2.41 CV 21.01 24.15 23.70 105.17 26.87 38.26 34.53 38.26 O-Ru Mean 541.61 155.90 347.27 38.42 338.33 449.02 98.75 7.05 Min 266.30 32.60 204.80 0.40 190.00 230.00 58.44 3.61 Max 698.50 213.20 480.20 147.90 559.00 715.00 153.60 11.23 SD 103.67 42.25 74.90 41.11 89.32 128.78 26.67 2.02 CV 19.14 27.10 21.57 107.01 26.40 28.68 27.01 28.68 AH Mean 654.04 153.31 464.58 36.13 426.00 291.65 79.65 4.58 Min 270.90 55.20 214.30 0.05 162.70 120.00 51.81 1.88 Max 763.50 274.60 591.90 130.70 660.60 779.00 167.40 12.23 SD 95.48 37.74 79.83 33.13 94.55 132.98 23.00 2.09 CV 14.60 24.62 17.18 91.70 22.20 45.59 28.88 45.59 CH Mean 680.96 177.63 455.31 47.99 474.82 276.15 68.27 4.34 Min 564.20 86.90 363.80 1.10 286.00 110.00 50.94 1.73 Max 779.20 218.20 577.00 157.90 635.00 460.00 86.81 7.22 SD 65.74 26.87 54.12 42.72 95.38 99.24 10.09 1.56 CV 9.65 15.13 11.89 89.02 20.09 35.94 14.78 35.94 Table A4. Characteristics of the tree layer and cover of shrubs, grasses and forbs in the stands selected for the study. STAND Tree Density (trees/ha) Tree Canopy Cover (%) Tree Main Species Shrubs * (%) Grasses (%) Forbs (%) 0-Bl (1) Boimorto 1000 90 Quercus robur L1: 60 L2: 10 30 20 (2) Cerceda 1000 80 Quercus robur, Castanea sativa L1: 30 L2: 25 35 10 (3) Cerqueiras 400 70 Quercus robur L1: 65 L2: 20 70 30 (4) Fragavella 500 70 Quercus robur L1: 25 L2: 40 50 40 (5) Monfero 500 70 Quercus robur L1: 20 L2: 30 60 20 O-My (6) A Lama 500 80 Quercus robur L1: 85 L2: 75 75 15 (7) Lovios 600 80 Quercus roburFagus sylvatica L1: 50 L2: 80 10 20 (8) Xestoso 900 70 Quercus robur L1: 75 L2: 60 10 5 O-Ru (9) A Rua 700 80 Quercus robur L1: 50 L2: 75 5 15 (10) Lourizán 200 70 Quercus robur L1: 30 L2: 60 35 25 AH (11) Arnuide —– —– —– L1: 25 L2: 70 40 —– (12) Cerqueiras 1100 60 Pinus pinaster L1: 25 L2: 60 60 20 (13) Cerponzóns 300 30 Pinus pinaster L1: 45 L2: 60 35 —– (14) Coto do Muiño 1200 40 Eucalyptus globulus L1: 60 L2: 75 25 35 (15) Fragavella —– —– —– L1: 70 L2: 50 35 20 (16) Oia 300 60 Pinus pinaster L1: 20 L2: 60 50 10 (17) Pouzos —– —– —– L1: 60 L2: 35 35 35 (18) Silleda —– —– —– L1: 70 L2: 35 65 —– CH (19) Eirexa —– —– —– L1: 30 L2: 50 25 10 (20) Rubiás —– —– —– L1: 30 L2: 60 35 15 * Shrubs are presented in two differentiated layers: L 1 : shrub stratum over 50 cm of height; L 2 : shrub stratum lower than 50 cm of height. Sustainability 2021,13, 13581 20 of 22 References 1. Silva-Pando, F.J.; Rozados Lorenzo, M.J.; González-Hernández, M.P. Grasslands and Scrublands in the Northwest of the Iberian Peninsula: Silvopastoral Systems and Nature Conservation. In Pasture Landscapes and Nature Conservation; Redecker, B., Finck, P., Härdtle, W., Riecken, U., Schröder, E., Eds.; Springer: New York, NY, USA; Berlin/Heidelberg, Germany, 2002; pp. 271–283. 2. Celaya, R.; García, R.R.; Benavides, R.; García, U.; Osoro, K.; Jáuregui, B.M.; García, R.R. Changes in heathland vegetation under goat grazing: Effects of breed and stocking rate. Appl. Veg. Sci. 2010,13, 125–134. [CrossRef] 3. Fuhlendorf, S.D.; Engle, D.M.; Elmore, R.D.; Limb, R.F.; Bidwell, T.G. Conservation of Pattern and Process: Developing an Alternative Paradigm of Rangeland Management. Rangel. Ecol. Manag. 2012,65, 579–589. [CrossRef] 4. González-Hernández, M.; Mouronte, V.; Romero, R.; Rigueiro-Rodríguez, A.; Mosquera-Losada, M.R.O. Plant diversity and botanical composition in an Atlantic heather-gorse dominated understory after horse grazing suspension: Comparison of a continuous and rotational management. Glob. Ecol. Conserv. 2020,23, e01134. [CrossRef] 5. López, C.L.; García, R.R.; Ferreira, L.; García, U.; Osoro, K.; Celaya, R. Impacts of horse grazing on botanical composition and diversity in different types of heathland. Rangel. J. 2017,39, 375–385. [CrossRef] 6. Rigueiro-Rodríguez, A.; Mouhbi, R.; Santiago-Freijanes, J.J.; Hernández, M.D.P.G.; Mosquera-Losada, M.R. Horse grazing systems: Understory biomass and plant biodiversity of a Pinus radiata stand. Sci. Agricola 2012,69, 38–46. [CrossRef] 7. González-Hernández, M.P.; Silva-Pando, F.J. Grazing effects of ungulates in a Galician oakwood (NW Spain). For. Ecol. Manag. 1996,88, 65–70. [CrossRef] 8. Castro, M. Silvopastoral Systems in Portugal: Current Status and Future Prospects. In Agroforestry in Europe; Rigueiro-Rodríguez, A., McAdam, J., Mosquera-Losada, M.R., Eds.; Springer: Dordrecht, Germany, 2008; Volume 6, pp. 111–126. [CrossRef] 9. Pantera, A.; Papadopoulos, A.; Papanastasis, V.P. Valonia oak agroforestry systems in Greece: An overview. Agrofor. Syst. 2018 , 92, 921–931. [CrossRef] 10. Pardini, A.; Nori, M. Agro-silvo-pastoral systems in Italy: Integration and diversification. Pastoralism 2011,1, 26. [CrossRef] 11. Varga, A.; Demeter, L.; Ulicsni, V.; Öllerer, K.; Biró, M.; Babai, D.; Molnár, Z. Prohibited, but still present: Local and traditional knowledge about the practice and impact of forest grazing by domestic livestock in Hungary. J. Ethnobiol. Ethnomed. 2020 ,16, 1–12. [CrossRef] 12. Castro, M.; Teixeira, A.; Fernández-Núñez, E. The nutritive value of different Mediterranean browse species used as animal feeds under oak silvopastoral systems in Northern Portugal. Agrofor. Syst. 2021,95, 269–278. [CrossRef] 13. Osoro, K.; Vassallo, L.M.; Celaya, R.; Martínez, A. Livestock production systems and the vegetation dynamics of Less Favoured Areas (LFAs): Developing viable systems to manage semi-natural vegetation in temperate LFAs in Spain. In Livestock Production in the European less Favoured Areas: Meeting Future Economic, Environmental and Policy Objectives through Integrated Research; Laker, J.P., Milne, J.A., Eds.; Macaulay Land Use Research Institute: Aberdeen, UK, 1999; pp. 133–143. 14. García, R.R.; Fraser, M.D.; Celaya, R.; Ferreira, L.M.M.; García, U.; Osoro, K. Grazing land management and biodiversity in the Atlantic European heathlands: A review. Agrofor. Syst. 2012,87, 19–43. [CrossRef] 15. Lasanta, T.; Khorchani, M.; Pérez-Cabello, F.; Errea, P.; Sáenz-Blanco, R.; Nadal-Romero, E. Clearing shrubland and extensive livestock farming: Active prevention to control wildfires in the Mediterranean mountains. J. Environ. Manag. 2018 ,227, 256–266. [CrossRef] [PubMed] 16. Gavazov, K.S.; Peringer, A.; Buttler, A.; Gillet, F.; Spiegelberger, T. Dynamics of Forage Production in Pasture-woodlands of the Swiss Jura Mountains under Projected Climate Change Scenarios. Ecol. Soc. 2013,18, 38. [CrossRef] 17. González-Hernández, M.P.; Silva-Pando, F.J. Nutritional Attributes of Understory Plants Known as Components of Deer Diets. J. Range Manag. 1999,52, 132. [CrossRef] 18. Jáuregui, B.M.; Celaya, R.; García, U.; Osoro, K. Vegetation dynamics in burnt heather-gorse shrublands under different grazing management with sheep and goats. Agrofor. Syst. 2007,70, 103–111. [CrossRef] 19. Osoro, K.; Ferreira, L.M.M.; García, U.; Jáuregui, B.M.; Rosa García, R.; Celaya, R. Diet selection and performance of sheep and goats grazing on different heathland vegetation types. Small Rumin. Res. 2013,109, 119–127. [CrossRef] 20. Osoro, K.; Mateos-Sanz, A.; Frutos, P.; García, U.; Ortega-Mora, L.M.; Ferreira, L.M.M.; Celaya, R.; Ferre, I. Anthelmintic and nutritional effects of heather supplementation on Cashmere goats grazing perennial ryegrass-white clover pastures. J. Anim. Sci. 2007,85, 861–870. [CrossRef] [PubMed] 21. Hiemstra, S.J.; de Haas, Y.; Mäki-Tanila, A.; Gandini, G. Local Cattle Breeds in Europe Development of Policies and Strategies for Self-Sustaining Breeds; Wageningen Academic Publishers: Wageningen, The Netherlands, 2010. 22. Heitschmidt, R.K.; Stuth, J.W. Grazing Management. An Ecological Perspective; Timber Press Inc.: Portland, Oregon, 1991. 23. Caughley, G.; Sinclair, A.R.E. Wildlife Ecology and Management; Blackwell Science, Inc.: Oxford, UK, 1994. 24. Holechek, J.L.; Pieper, R.D.; Herbel, C.H. Range Management: Principles and Practices, 5th ed.; Prentice Hall: Upper Saddle River, NJ, USA, 2004. 25. Yang, Y.; Li, Q.; Garber, P.A.; Grueter, C.C.; Ren, G.; Wang, X.; Huang, Z.; Xiang, Z.; Xiao, W.; Behie, A. Cafeteria-style feeding trials provide new insights into the diet and nutritional strategies of the black snub-nosed monkey ( Rhinopithecus strykeri ): Implications for conservation. Am. J. Primatol. 2020,82, e23108. [CrossRef] [PubMed] 26. Tassone, S.; Fortina, R.; Valle, E.; Cavallarin, L.; Raspa, F.; Boggero, S.; Bergero, D.; Giammarino, M.; Renna, M. Comparison of In Vivo and In Vitro Digestibility in Donkeys. Animals 2020,10, 2100. [CrossRef] [PubMed] 27. Weiss, W.P.; Tebbe, A.W. Estimating digestible energy values of feeds and diets and integrating those values into net energy systems. Transl. Anim. Sci. 2018,3, 953–961. [CrossRef] [PubMed] Sustainability 2021,13, 13581 21 of 22 28. AFRC. Energy and protein requirements of ruminants. An Advisory Manual Prepared by the AFRC Technical Committee on Responses to Nutrients; CAB International: Wallingford, UK, 1993. 29. Maizeret, C.; Sung, T.M. Etude du régime alimentaire et recherche du déterminisme fonctionnel de la sélectivitéchez le chevreuil (Capreolus capreolus) des Landes de Gascogne. Gibier Faune Sauvag. 1984,3, 63–103. 30. Sineiro, F.; Osoro, K.; Díaz, N. Bases para la producción e intensificación ganadera en el monte gallego: La utilización de la vegetación espontánea y la siembra y mejora del pasto. In Pastos y Forrajes en Alimentación Animal; Tella, D.C., Ed.; Actas de la XXII Reunión Científica de la S.I.N.A.: Santiago de Compostela, Spain, 1984; pp. 195–219. 31. Putman, R.J.; Pratt, R.M.; Ekins, J.R.; Edwards, P.J. Food and Feeding Behaviour of Cattle and Ponies in the New Forest, Hampshire. J. Appl. Ecol. 1987,24, 369. [CrossRef] 32. Alibes, X.; Tisserand, J.L. Tableaux de la Valeur Alimentaire pour les Ruminants des Fourrages et Sous-Produits D’origine Méditerranéenne; CIHEAM: Zaragoza, Spain, 1990. 33. Yan, T.; Agnew, R.E. Prediction of nutritive values in grass silages: I. Nutrient digestibility and energy concentrations using nutrient compositions and fermentation characteristics. J. Anim. Sci. 2004,82, 1367–1379. [CrossRef] 34. Wan, H.F.; Chen, W.; Qi, Z.L.; Peng, P.; Peng, J. Prediction of true metabolizable energy from chemical composition of wheat milling by-products for ducks. Poult. Sci. 2009,88, 92–97. [CrossRef] 35. Stergiadis, S.; Allen, M.; Chen, X.; Wills, D.; Yan, T. Prediction of nutrient digestibility and energy concentrations in fresh grass using nutrient composition. J. Dairy Sci. 2015,98, 3257–3273. [CrossRef] [PubMed] 36. Rivas-Martínez, S.; Penas, Á.; González, T.E.D.; Cantó, P.; Del Río, S.; Costa, J.C.; Herrero, L.; Molero, J. Biogeographic Units of the Iberian Peninsula and Baelaric Islands to District Level. A Concise Synopsis. In The Vegetation of the Iberian Peninsula; Springer: Berlin/Heidelberg, Germany, 2017; pp. 131–188. [CrossRef] 37. Izco, J.; Amigo, J.; García-San León, D. Análisis y clasificación de la vegetación leñosa de Galicia (España). Lazaroa 1999 ,20, 29–47. 38. Goering, H.K.; Van Soest, P.J. Forage fiber analysis (aparatus, reagents, procedures and some applications). Agr. Handb. Agr. Res. Serv. USDA 1970,379, 1–19. 39. Tilley, J.M.A.; Terry, R.A. A TWO-STAGE TECHNIQUE FOR THE IN VITRO DIGESTION OF FORAGE CROPS. Grass Forage Sci. 1963,18, 104–111. [CrossRef] 40. Alexander, R.H. The Establishment of a Laboratory Procedure for the " in vitro " Determination of Digestibility; Research Bulletin n ◦ 42; The West of Scotland Agricultural College: Perth, UK, 1969. 41. SAS Institute. SAS/STAT®9.1 User’s Guide; SAS Institute, Inc.: Cary, NC, USA, 2004. 42. Draper, N.R.; Smith, H. Applied Regression Analysis, 3rd ed.; John Wiley & Sons: New York, NY, USA, 1998. 43. Neter, J.; Wasserman, W.; Kutner, M.H. Applied Linear Statistical Models: Regression, Analysis of Variance and Experimental Designs; 3rd ed.; McGraw-Hill: Boston, MA, USA, 1990. 44. Kozak, A. Effects of multicollinearity and autocorrelation on the variable-exponent taper functions. Can. J. For. Res. 1997 ,27, 619–629. [CrossRef] 45. Schröder, J.; Rodríguez, R.; Vega, G. An age-independent basal area increment model for maritime pine trees in northwestern Spain. For. Ecol. Manag. 2002,157, 55–64. [CrossRef] 46. Van Soest, P.J.; Jones, L.H.P. Effect of silica in forages upon digestibility. J. Dairy Sci. 1968,51, 1644–1648. [CrossRef] 47. Smith, G.S.; Nelson, A.B.; Boggino, E.J. Digestibility of Forages In Vitro as Affected by Content of “Silica”. J. Anim. Sci. 1971 ,33, 466–471. [CrossRef] [PubMed] 48. Minson, D.J. Effect of chemical composition on feed digestibility and metabolizable energy. Nutr. Abstr. Rev. 1982,52, 591–615. 49. Varhegyi, I.; Szentmihalyi, S.; Varhegyi, J.; Simon, Z. Relation of crude fibre content and cell-wall constituents to dry matter digestibility in roughages. Acta Agron. Hung. 1987,36, 341–350. 50. Paton, D. Elaboration of a multi-variate model for the determination of the metabolizable energy of Mediterranean bushes based on chemical parameters. J. Arid. Environ. 2003,53, 271–280. [CrossRef] 51. Mahipala, M.K.; Krebs, G.; McCafferty, P.; Dods, K.; Suriyagoda, B. Faecal indices predict organic matter digestibility, short chain fatty acid production and metabolizable energy content of browse-containing sheep diets. Anim. Feed. Sci. Technol. 2009 ,154, 68–75. [CrossRef] 52. Nagy, J.G.; Hakonson, T.; Knox, K.L. Effect of quality on food intake in deer. Trans. of the Thirty-Fourth, N. Amer. Wildl. and Nat. Res. Conference; Wildlife Management Institute: Washington, DC, USA, 1969; pp. 146–154. 53. Dro˙ zd˙ z, A.; Osiecki, A. Intake and digestibility of natural feeds by roe-deer. Acta Theriol. 1973,18, 81–91. [CrossRef] 54. Crawford, H.S. Seasonal Food Selection and Digestibility by Tame White-Tailed Deer in Central Maine. J. Wildl. Manag. 1982 ,46, 974. [CrossRef] 55. Van Soest, P.J. Rice straw, the role of silica and treatments to improve quality. Anim. Feed. Sci. Technol. 2006 ,130, 137–171. [CrossRef] 56. Goddard, J.; McLean, E. Acid-insoluble ash as an inert reference material for digestibility studies in tilapia, Oreochromis aureus. Aquaculture 2001,194, 93–98. [CrossRef] 57. McGeough, E.J.; O’Kiely1, P.; Kenny, D.A. A note on the evaluation of the acid-insoluble ash technique as a method for determining apparent diet digestibility in beef cattle. Ir. J. Agric. Food Res. 2010,49, 159–164. 58. Sales, J.; Janssens, G.P.J. Acid-insoluble ash as a marker in digestibility studies: A review. J. Anim. Feed Sci. 2003 ,12, 383–401. [CrossRef] Sustainability 2021,13, 13581 22 of 22 59. Heinrichs, S.; Schmidt, W. Dynamics of Hedera helix L. in Central European beech forests on limestone: Results from long-term monitoring and experimental studies. Plant Ecol. 2014,216, 1–15. [CrossRef] 60. Van Uytvanck, J. The Role of Large Herbivores in Woodland Regeneration Patterns, Mechanisms and Processes. Ph.D Thesis, Research Institute for Nature and Forest, Ghent University, Ghent, Belgium, 2009.