Estimating Energy Concentrations in Wooded Pastures of NW Spain Using Empirical Models That Relate Observed Metabolizable Energy to Measured Nutritional Attributes

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² = 0.64) is applied, as opposed to the use of the specific prediction equations for each vegetation type (R² = 0.61, 0.66, 0.71 for oak woodlands; R² = 0.70 heather-gorse dominated heathlands, R² = 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.

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,2,3,4]. Horse grazing, for example, can decrease gorse dominance and promote heath composition in Atlantic heathlands [4,5,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,8,9,10,11,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 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,23,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,30,31,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,34,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.

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′ N, 7°51′ 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 (10 sites) and two types of heathlands: Atlantic heath-gorse 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

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.

and

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).

(Appendix A).

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 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=\frac{1}{1-R_i^2}$$

where R²_i is the multiple determination coefficient obtained when the ith independent variable X_i is 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 ($\overline{E}$), mean absolute error (|$\overline{E}$|), root mean square error (RMSE) and the adjusted coefficient of determination (R²_adj). These expressions can by summarized as follows:

$$Bias=\overline{E}={\displaystyle \sum }_{i=1}^n\frac{\left(y_i-\widehat{y_i}\right)}{n}$$

$$MAE=\left|\overline{E}\right|={\displaystyle \sum }_{i=1}^n\frac{\left|y_i-\widehat{y_i}\right|}{n}$$

$$RMSE=\sqrt{\frac{{\sum }_{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2}{n-1} }$$

$$R_{adj}^2=1-\frac{{\sum }_{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2·\left(n-1\right)}{{\sum }_{i=1}^n{\left(y_i-\overline{y}\right)}^2·\left(n-p\right)}$$

where $e_i=y_i-{\stackrel{⌢}{y}}_i$ and y_i, ${\widehat{y}}_i$ and $\overline{y}$ are the measured, predicted and average values of the dependent variables, respectively; n is the total number of observations used to fit the model and p is 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² 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, 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 Table A1, Table A2,

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

and

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₁: shrub stratum over 50 cm of height; L₂: shrub stratum lower than 50 cm of height.

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.

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. 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 1	Equation Code	Equations	R2	MEF	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. ¹ Pooled data = all data for the 5 vegetation types

. 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 adding lignin and silica (as well as their natural logarithm transformations) as predictor variables, and final R² reached values between 0.61–0.71. In O-Ru and O-My, the final R² 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² 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. 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	R2	MEF	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.

). 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² values for Summer (0.57) and better adjustments for Spring and Winter (R² 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. 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	R2	RMSE%	% 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 type
	O-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

). 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 (Figure 3, Figure 4 and Figure 5).

A clear tendency towards overestimation was observed for the plant group of heaths with the use of the general model (Figure 4A) that is partially corrected with the use of the specific models (Figure 4B) for each vegetation type (i.e., Equations (1c), (2e), (3c), (4d) and (5b)). This same tendency can be followed by season in Figure 4C.

4. Discussion

4.1. Establishment of the Prediction Equations

In the current study, we found that the content of ADF in combination with lignin or silica, as well as with DM, can be optimal predictors of ME for the understory of oak woodlands; in the case of heathlands, CP joined ADF and silica as a significant ME predictor. This is in agreement with a variety of research conducted on the effect of chemical composition on forage digestibility and metabolizable energy, which showed that fiber content and CP can be optimal predictors of ME [35,46,47,48,49,50,51]. Because digestibility and CP content are important components of quality in diets, whereas fiber is usually detrimentally associated with digestibility [52,53,54], CP is often positively associated with ME while fiber is negatively related with ME in the equations [35,46,47,50]. Phuong et al., (2013) reported that fibrous parameters of feed, especially ADF (but also NDF), reduce energy-use efficiency because they contain limited energy to sustain microbial growth, decrease availability of nutrients in high fiber diets and tend to accelerate rumen filling and limit feed intake [35]. The relationships of ADF, silica and CP with ME in the prediction equations here obtained showed a consistency with those findings. Nevertheless, the positive relation of lignin with ME in our study was not coherent with the fact of it being an undigestible fraction in plants. According to Van Soest (2006), correlation of lignin contents and degradability manifest poorly across a large amount of research papers, and it is not uncommon that the association of lignin with digestibility can be difficult to demonstrate [55].

Studies on the prediction of the digestibility of forages (and therefore extended to the prediction of forage energy concentrations) underline stepwise regression analyses in which the digestibility of organic matter is most accurately predicted by equations that include cell-wall constituents and acid-detergent lignin and silica [35,46,47]. The relationship between silica and the digestibility of forages showed that silica can be a major factor in limiting nutritional value [46], and the compilation of studies revealed an evident negative association of silica with IVOMD [55]. An average decline of 3.0 units of digestibility per unit of silica in the dry matter was found in grasses [46], and Smith et al. (1971) reported that silica decreases the digestibility of organic matter by approximately 1% for each percentage increase in the silica content of the forage [47]. This negative relationship between silica and digestibility led to the consideration of acid-insoluble ash (mainly silica) as a good internal marker to estimate dietary digestibility in mammals and poultry (Van Keulen and Young, 1977, Vogtmann et al., 1975) [56]. The use of other internal markers, such as crude fiber and hydrolysis-resistant organic matter, have also been recommended for measurement of apparent digestibility coefficients (De Silva and Anderson, 1995) [56]. Our results showed a noticeable increase of 0.20–0.30 units in R² for oak woodland prediction equations (Equations (2b) and (3b)) when DM was entered as a predictor in the equations. Silica was incorporated later in O-My and in Atlantic heathlands, improving slightly the accuracy of the ME final equations. Although the contribution of silica in the raising of R² in our prediction equations was small, this parameter should not be overlooked, because estimates of dietary digestibility obtained using the acid-insoluble ash method were in tolerable agreement with those obtained using the standard total fecal collection method [57], which represents a laborious method of determination of in vivo diet digestibility.

The success of acid-insoluble ash as a suitable digestibility marker has depended on the amount of silica content in the forage, with recommendations of a minimum of 7.5 g/kg on a dry matter (DM) basis to get accurate measurements [58]. Based on our silica results, most of the plant species included in the prediction equations surpass that mean value, except some forbs such as Asphodelus albus (associated with forest stands previously burnt) and grasses like Agrostis curtisii. Therefore, if those ME equations with silica as a predictor were to be applied to vegetation types where these plant species are very dominant, some caution should be taken, and perhaps alternative equations with no silica as a predictor could be more suitable (i.e., Equations (4b) and (2c)). Finally, the inclusion of silica in the equations could have some interesting management implications as acid-insoluble ash has been also used to determine feed intake in grazing animals [58], a factor that is related to the continuance or permanence of pasture.

4.2. Accuracy and Applicability of the Prediction Equations

In predicting a numerical outcome with a statistical model, the residuals provide an idea of the differences between predicted and observed values (i.e., in our study values of ME using the AFRC formula with IVOMD as a sole predictor). The analysis of the residuals to test for the accuracy of the prediction equations obtained in this study revealed that the use of the general model (Equation (6d) from pooled data) as opposed to the use of the specific prediction equations for each vegetation type (i.e., Equations (1c), (2e), (3c), (4d) and (5c)) or for each season, would mean a loss of precision in the ME estimation (Table 3). Accuracy, expressed as a combination of RMSE and bias and MAE values (a measure of errors between paired observations), was higher when the general model was applied, and the magnitude of the loss of precision was superior for Atlantic heathlands (RMSE increase of 29%). ME estimation in continental heathlands and oak woodlands would lose 1 and 11% of accuracy, respectively, if the general model were chosen. On the other hand, the degree of total loss of accuracy in predicting seasonal ME with the general model as opposed to the use of the specific equation for each season would be much smaller (5.4%). Only Winter ME estimates are expected to be substantially influenced by the use of one or another (the expected loss of precision being around 25%). RMSE increased by only 1.8% for Spring prediction when applying the general model, whereas for Fall and Summer the precision loss with the general model is expected to be around 4%. Because the quality of the fit does not necessarily reflect the quality of the prediction, the cross-validation of each model was examined. The equivalent estimate to the R² of the fitting phase, namely MEF in the cross-validation, followed a similar tendency, and showed again for those models fitted specifically for each vegetation type higher values than those where the general model was applied. Therefore, the quality of the prediction is expected to improve when using the specific prediction equations as opposed to the general model. Conversely, however, cross-validation performed for each season revealed MFE values approaching closely those obtained from applying the general model (and only still differentiated for Winter), thus confirming once again that only Winter ME estimates are expected to be substantially improved by the use of the specific season model as opposed to the general model. Shrubs, and especially those in the Ericaceae family, show a nutritional composition characterized by a very small seasonal variation [17].

The range of predicted ME concentrations in oak woodlands and heathlands narrowed with respect to the observed ME concentrations because the minimum values of the range increased (+0.6 for O-Ru, +1.2 for O-Bl, +0.8 for both heathlands) and maximum values were still similar to observed ME (around 13 MJ/kg DM). The results of graphing the observed ME against the predicted ME revealed the degree of dispersion of the data when using a prediction equation based on nutritional parameters (ADF, lignin, silica, DM, and CP) as opposed to the one based solely on IVOMD. Our attempt to determine whether plant groups (forbs, grasses, heaths, other shrubs and trees) were differently affected by the usual discrepancy between observed ME and predicted ME showed a clear tendency towards the overestimation of heaths when using the general model. Such overestimation could be partially corrected with the use of the specific models for each vegetation type. This is important to take into account in the case of ME predictions performed with the general model for heathlands that are highly dominated by heaths, because the observed ME data obtained from AFRC (1993) could already be overestimating the actual ME concentration of heaths. Yan and Agnew (2004) reported an over-prediction of ME concentration of 1.2 MJ/kg DM when using AFRC ME estimations for poor quality silages, the under-prediction for very good quality silage being 1.3 MJ/kg DM. As an example of possible implications for management, efficient reducers of fuel load like goats have ME requirements for maintenance around 7.4 MJ/kg DM (0.424 MJ/kg BW ^0.75, taken from NRC 1981), for an average animal of 45 kg. In the case of beef cattle, similar values of 7.6 MJ/kg DM have been approached for Winter season, raising in general as an average of 8.6 MJ/kg DM (NRC 1996) for a mature animal of 1000 lb (453,592 kg). In both continental and Atlantic heathlands, the observed ME average values were under the minimum requirements for maintenance because of a high percentage of heaths and Halimium spp., and such overestimation could result in a substantial divergence facing management decisions (for example of the optimum animal stoking rates).

There is little literature discussing ME concentrations of browse considered in our study. Most data on pastures’ ME concentrations are relative to fresh grass forages with ME values between 7.5–12.6 MJ/kg DM [28,35], or grass silages with ME around 13.6 MJ/kg DM [33]. Over 2/3 of the plants included in Atlantic and continental heathlands showed observed ME values below 5 MJ/kg DM, which are similar observed ME figures to those observed in the Ericaceae plants reported by Paton (2003) [50]. Other groups of plants studied here showed similar ME concentrations to those found in the same study. For example, gorse such as U. minor and Ulex europaeus included in Atlantic heathlands reached ME concentrations of 8.5 MJ/kg DM and between 6.9–7.2 MJ/kg DM, respectively. Gorse is highly selected by horses in heathlands, decreasing its biomass and dominance and leading to the reduction of fire risk as well as increasing plant diversity [4,5,6,31]. Unlike in heathlands, observed ME concentrations in the three types of oak woodlands revealed more records with values superior to 7 MJ/kg DM. ME concentrations of around 8 MJ/kg DM obtained for the leguminous shrubs Cytisus striatus and Genista florida in our study were similar to those found in plants of the same genera in Castro et al. (2021) [12]. The highest number of plants reaching minimum ME concentration for maintenance occurred in O-Bl, and were mainly associated with the shrubs Hedera helix, Lonicera periclymenum and Ilex aquifolium; to the forbs Asphodelus albus and Teucrium scorodonia and to grasses. These are relevant plant species highly preferent in roe deer diets [59,60], and observed ME concentrations in the present study are also in the range reported by Paton (2003) [50]. Some forbs in O-Bl and O-My showed higher ME concentrations because of superior nutrient content and better digestibility, and were close to reaching in Spring the minimum ME requirements of goats for late pregnancy (13 MJ/kg DM). This would balance the low ME values provided by other wooded plant species.

The stepwise regression analysis of the present study generated prediction equations of ME based on ordinary nutritional parameters that are easily obtained from laboratories (i.e., ADF, lignin, crude protein, dry matter), which could facilitate ME estimations in grazing management scenarios where absolute accuracy is not essential. This could greatly simplify grazing management decisions in an extensive forest grazing setting with a wide variety of plants, being an alternative to in vivo digestibility trials and in vitro digestibility analysis, which are not always available and are more difficult to conduct.

There has been a growing interest in predicting the energy value of diets by establishing relationships with chemical parameters that are easily analyzed in laboratories, but, to the best of our knowledge, this is the first time that regression equations are presented in order to predict the ME of the main Atlantic pedunculate oak woodlands and heathlands.

The wide-ranging information on the nutritional characteristics of the plant species from which the ME predictive equations have been constructed provides a reference source when these equations are to be applied in similar vegetation types. Prediction equations in the current study have been developed by using a great range of browse quality, also including some seasonal variation of nutritional characteristics, so that they can be broadly applicable in similar vegetation types depending on the available data on chemical composition. Our results showed that the prediction accuracy of ME is reduced when the general model is applied, as opposed to the specific prediction equations for each vegetation type, especially in the case of Atlantic heathlands. The general model 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.