Mean species cover: a harmonized indicator of shrub cover for forest inventories

Because shrub cover is related to many forest ecosystem functions, it is one of the most relevant variables for describing these communities. Nevertheless, a harmonized indicator of shrub cover for large-scale reporting is lacking. The aims of the study were threefold: to define a shrub indicator that can be used by European countries for harmonized shrub cover estimation using data from their respective national forest inventories (NFIs); to quantify the effects of using different NFI field cover scales; and to establish bridges to facilitate harmonized estimation. Data for shrub species cover from the Third Spanish NFI together with scales for cover assessment from 16 European NFIs were used. The indicator, mean species cover (MSC), was defined for each species and each European forest category. Estimates of MSC calculated using species covers recorded for field plots, with 1% interval widths (MSCobs), were compared with the MSC values that would be obtained for the same data with the different European cover scales (MSCpred). Residuals calculated as differences between MSCobs and MSCpred were analyzed, and a linear mixed model was used as bridging function to adjust predictions and thus further harmonize estimates. Scales with only two or three intervals produced the greatest residuals, while all the other analyzed scales had residuals less than 5%. Most scales, except those most similar to Braun-Blanquet, displayed a tendency to be unreliable for larger covers. The proposed mean species cover indicator provides comparable estimates for shrub communities at large scales. The linear models improved the harmonization of MSC for the scales having two and three intervals.


Introduction
Analyzing non-tree vegetation is important for forest management, whether it is aimed at conservation or production (Hart and Chen 2006;Nilsson and Wardle 2005). Understory vegetation, and shrubs in particular, provide shelter and food for fauna (Carrilho et al. 2017;Fortuny et al. 2014;Mangas et al. 2008), form communities of great ecological value such as in floodplains, act as bioindicators for erosion (Andreu et al. 1998;Francis and Thornes 1990) and fire (Cerdà and Doerr 2005;Fréjaville et al. 2016), and contribute to the recovery of deforested and degraded landscapes (Tasser and Tappeiner 2002;Valdecantos et al. 2009). Nontree forest vegetation is also important as a living carbon sink (Peñuelas et al. 2002) that has been increasingly recognized by the Intergovernmental Panel on Climate Change (IPCC) (IPCC 2014).
Shrub cover, along with species composition, is among the most relevant attributes for describing shrub communities. Shrub cover reflects the availability of soil, water, and nutrients that a plant can use to create biomass. It influences infiltration and potential erosion, wildlife habitat, and the availability of forage (Muir and McClaran 1997). Also, the degree to which shrub cover influences regeneration and survival increases or decreases depends on whether speciesspecific positive and negative effects are in balance (Nilsson and Wardle 2005;Padilla and Pugnaire 2007). Specifically, while shrub cover positively affects tree seedling survival and initial growth by contributing to water and temperature balances (Gómez-Aparicio et al. 2004), it can also negatively affect seedling growth via competition or interference resulting from sharing limited resources or release of chemicals that harm nearby plants (Kitzberger et al. 2000;Padilla and Pugnaire 2006).
In addition to their local ecological importance, the importance of shrubs has been recognized by multiple international reporting organizations. First, as parties to both the Kyoto Protocol and the United Nations Framework Convention on Climate Change, European countries individually and in aggregate as the European Union (EU) report greenhouse gas emissions for the major Land Use, Land Use Change and Forestry carbon pools of which one, live biomass, includes shrubs. Second, for the 2015 Global Forest Resources Assessment (FRA), FAO has defined other wooded land (OWL) to include land whose "combined cover of shrubs, bushes and trees" is greater than 10% (FAO 2012a). The 2015 FRA additionally reports estimates of area and growing stock for OWL for all European countries (FAO 2015). Third, Forest Europe (previously Ministerial Convention on the Protection of Forests in Europe, (MCPFE 1998)) uses the same OWL definition as FAO and includes the area of this land cover category as Indicator 1.1, Forest area and OWL, for Criterion C1, Maintenance and appropriate enhancement of forest resources and their contribution to global carbon cycles (Forest Europe 2015).
For reporting at cross-regional scales, national forest inventories (NFI) are among the most important sources of forest and forest-related information because of their large numbers of variables and large numbers of field plots (Tomppo et al. 2010). Apart from trees, shrubs are the vegetation component for which the most complete and detailed information is available from NFIs. Chirici et al. (2012) reported that within the framework of NFIs, species composition and cover are the most appropriate indicators for characterizing shrubs.
For NFI purposes, Vidal et al. (2016a) defined shrub cover as the "horizontal projection of the outermost crown limits of a shrub or a group of shrubs." When applied to a single species, the result it is called "species cover" and when applied without distinctions among species is called "total cover" (Wilson 2011). Although shrub cover could be assessed by the area it occupies, it is generally assessed as an area percentage in the range 0-100%. For estimating shrub cover, NFIs use a variety of cover scales (Bonham 1989), (Tables 1, 2), although most countries use a modified Braun-Blanquet scale (Braun-Blanquet 1965) based on visual sampling, mostly without consideration of density (Wilson 2011). Moreover, not all shrub species are monitored by the NFIs of all countries; rather, the majority of NFIs use species lists (Alberdi et al. 2010) related to their inventory aims and the understory structure and composition of their forests . Thus, consistent, crosscountry reporting is extremely difficult. The 2015 Global FRA, FAO addressed this issue, not only for shrubs but all variables, by directing that 'if the national categories differ substantially from the FRA categories, countries should try reclassifying the national data to the FRA categories" (FAO 2012b). However, despite the importance of shrubs for both ecological and international assessment purposes, and despite the FRA emphasis, a harmonized indicator for large-scale shrub assessment is lacking.
Because the natural environment is structured in both space and time, and organisms respond to this structure (McGarigal et al. 2016), any shrub assessment should reflect this structure (Bastos et al. 2016). The European Forest Categories (EFC) (EEA 2006) were introduced to improve the quality of information reported by quantitative indicators of the State of Europe's Forests 2011 (Forest Europe, UNECE and FAO 2011). They are the most applicable forest-type classification for NFIs nowadays , and, therefore, their influence on any harmonized shrub estimation should be considered (Table 3).
The aims of the study were threefold: (1) to define a shrub indicator that can be used by European countries for harmonized shrub cover estimation using data from their respective NFIs, (2) to quantify the effects of using different NFIs field cover scales on harmonized estimates, and (3) to establish bridges to facilitate harmonized estimation.

Data
The dataset used for this study consisted of observations for 64,221 sample plots from the Third Spanish national forest inventory (SNFI-3). These sample plots were distributed throughout the Iberian Peninsula (excluding Portugal) and were classified with respect to the 14 EFCs. Table 3 shows the number of species in each EFC.
The SNFI has established permanent sample plots at the nodes of a 1-km × 1-km grid. SNFI-3 was conducted between 1997 and 2007 and included shrub measurements for 83,010 sample plots of which 18,789 plots were not considered for this study because of EFC classification difficulties, mostly for mixed stands. Shrub attributes were measured for circular, 10-m radius, sample plots. The SNFI visually assesses mean height (h) and cover for each shrub species using a percentage scale with 1% interval widths.
The Spanish NFI shrub assessment is based on shrub taxa lists defined using criteria based mainly on shrub dominance in the defined NFI forest stratum. However, there are also non-dominant species selected as bioindicators or key species. We deal with data from 190 different shrub taxa (mainly species) that could occur in one or more of the EFCs. For this study, we considered only the well-represented shrub species, i.e., those recorded for at least five sample plots in an EFC; a total of 152 species was considered, being 822 the sum of the of all possible shrub species appearing in each EFC (Table 3). We removed 38 species because most European countries include only dominant species in their shrub taxa lists ).

Definition of Mean Species Cover Indicator
With the aim of defining an indicator that could be used by most European countries, the different attributes measured by European NFIs were considered (Alberdi et al. 2010). Shrub cover was the non-tree variable most commonly recorded by European NFIs, and therefore indicators based on this variable have the greatest possibility for harmonization.
Moreover, a preliminary study of the SNFI-3 dataset analyzed the particular features of shrub cover for different species and the effect of NFI scales on plot shrub cover assessments. As an example, we include analyses for three species that are characteristic of forest types or climatic conditions that are well represented in Spain: Calluna vulgaris L., Cistus ladanifer L. and Erica arborea L.
Taking into account the cover distribution differences among species and EFCs, we defined the indicator mean shrub species cover by European Forest Category (MSC), as the mean shrub cover for plots where the target species was present. The number of plots where the species was not present should also be mentioned as important complementary information, thereby providing observations for two values, occurrence and MSC value. NFIs that collect sufficient field measurements to estimate the proposed indicator, together with their corresponding cover scales, are shown in Table 2.

Mean species cover estimation and scale effect analysis
For each species in each EFC, MSC was estimated using observed field data recorded during SNFI-3 using a percentage scale denoted S0 with 1% interval widths. This MSC estimate is hereafter characterized as "observed MSC" (MSCobs).
MSC was also estimated for the sample plots using each of the different European NFIs scales (Table 1). These estimates are characterized as "predicted MSC" (MSCpred). For a specific cover scale, MSCpred was calculated by assigning each plot to the midpoint of the interval that included the Table 1 NFI cover scales for shrubs used by the countries "r" denotes rare species, and "+" denotes scarce species Norway, Finland and Lithuania are not measuring shrub cover in the NFI survey observation of the plot's species cover (Faber-Langendoen et al. 2007;Canullo et al. 2011).
Because the predicted values were obtained from the observed data (in a percentage scale with 1% interval widths, S0), species with covers of less than 1% were assigned a cover of 1%, assuming that all NFIs are able to distinguish between the presence and absence of each species.
For each species in each EFC, the indicator MSCpred was estimated for the six cover scales (S1, S2, S3, S4, S5 and S6, Table 1) and compared with the MSCobs. The residuals calculated as differences between MSCobs and MSCpred were calculated and analyzed. In addition, mean error (ME, Eq. 1), mean absolute error (MAE, Eq. 2) and root mean square error (RMSE, Eq. 3) were also calculated for each scale: where MSCobs sp and MSCpred sp were the observed and predicted MSC for the species sp, and n was the number of species in each EFC.

Bridging function
Given the results from the previous section and with the objective of improving MSC harmonization, a bridging function based on a model of the relationship between MSCobs and MSCpred was developed. The model was fit for each cover scale using the MSC values of the species observed for at least five field plots (Table 3). Considering that differences in the accuracy of predictions could be related to the EFC, a mixed model with random effects corresponding to EFCs was used. The use of mixed models would allow the application of bridging functions for correcting MSC whether plot is located in a known EFC or not. For each scale, we formulated a common linear mixed model, where MSCobs sp j was the observed MSC for the species sp in EFC j, MSCpred sp j was the predicted MSC for the same species and EFC, and a i and a ij were, respectively, the fixed and random coefficients to be estimated.
With the aim of determining whether the bridging function could be applied to an independent dataset, we randomly selected a subsample of half of the data, i.e., 411 pairs of observed and predicted MSC sp j as calibration data for fitting the model and the remaining 411 as validation data for assessing accuracy. The model was fit to the calibration data and then applied to the validation data using MSCpred to obtain "corrected MSC" (MSCcorr). Residuals between MSCobs and MSCcorr together with ME, MAE and RMSE using Eqs. 1, 2 and 3, were calculated. The process of splitting the data into calibration and validation datasets, fitting the models, and calculating the accuracy measures was repeated 100 times, and appropriate summary statistics were calculated.
Using all the data, final models were fit using the restricted maximum likelihood method (REML). Predictor variables were assessed to determine whether they contributed to statistically significantly improving the quality of fit of the model to the data at the p = 0.05 level. Goodness-of-fit measures for the linear mixed models were also calculated. For this study, the lmmR2 and lme procedures in R (2014) were used, although other software could also be used.

Results
Firstly, it is interesting to note that the empirical cumulative distribution functions of covers recorded in field plots were different among species (or taxon) and even among EFCs for the same species (Fig. 1). Two extreme situations were found: (4) MSCobs sp j = a 0 + a 0j + a 1 + a 1j ⋅ MSCpred sp j + sp j , • Some species had smaller cover values for most of the plots of an EFC; for instance, Cistus ladanifer in EFC 8 ( Fig. 1-top) and Erica arborea in EFC 14 ( Fig. 1-bottom), resulting in distributions with steep slopes for the smaller covers. • The numbers of plots were approximately proportionally distributed over the entire cover range; for example, Cistus ladanifer in EFC 14 ( Fig. 1-top). For these species covers, the distributions were fairly linear.
These differences influence the resulting values of the MSCobs indicator, highlighting the importance of providing estimates separately by species and EFC (e.g., MSCobs of Cistus ladanifer in EFC 8 is 24%, while for EFC 14 is 36%). Table 3 summarizes the statistics for MSCobs in the  (Table 3). Mean species cover and occurrence, in brackets, of Cistus ladanifer are 24% (8%) and 36% (16%) for EFC 8 and EFC14, respectively, while for Erica arborea are 24% (20%), 27% (8%) and 16% (10%) for EFC 8, EFC 10 and EFC14, respectively different EFCs. Generally, the smallest MSC estimates were obtained for EFC 10 (8.6%), while the greatest estimates were obtained for EFC 13 (21.7%). Figure 2 shows how shapes of the cover cumulative distributions together with the NFI scales can influence the MSC estimates. As expected, scales with larger numbers of intervals in the smaller cover classes more accurately mimicked the distribution curves described in the previous section, whereas scales with smaller numbers of intervals such as S5 or S6 resulted in greater mean errors.

Scale effect on the mean species cover estimation
The residuals in the MSC estimates, i.e., MSCobs minus MSCpred, using different scales are shown in Fig. 3a. The scale with cover intervals of 10% width (S1) tended to underestimate the indicator (ME = 0.5), especially for species with large covers. The same effect, but less marked, can be seen for scale S2. However, the use of these scales produced no residuals greater than 5%.
The estimates obtained with scales S3 (modified Braun-Blanquet scale) and S4 were quite accurate with ME ≈ 0, and independent of the EFC (Fig. 4). The results shown in Fig. 4 for the other scales highlight the importance of estimating cover separately by EFC.
Scales S5 and S6, with only four and two cover intervals, respectively, overestimated MSC with the overestimation, especially pronounced for S6 where the ME is larger than 15%. However, while most of the residuals were less than 5% and MAE was 1.5 for scale S5, for scale S6 the residuals could be greater than 20% and MAE larger than 15%. The smallest differences were found for EFC13 for which the values for MSCobs were largest (Table 3).

Bridging function
The results of the iterative analysis for which the regressionbased correction was estimated for training datasets and applied to independent validation datasets showed that the mean errors always decreased after applying the bridging function (Table 4). Therefore, when comparing residuals calculated for MSCpred with those calculated for MSCcorr, ME not only approached zero, but clear reductions in MAE and RMSE were also achieved for all analyzed scales (Table 4). It is important to highlight that even the RMSE values for the most unfavorable scale, i.e., S6, were on average less than 5%.  Table 1). Observed MSC for Cistus ladanifer and Calluna vulgaris were, respectively, 26 and 14%. The residuals for estimates with scale S2 were less than 1% for both species, while for scale S6 residuals were, respectively, 5 and 7% Table 5 includes the final set of coefficient estimates together with marginal and conditional R 2 for each scale obtained when fitting regressions with all species and EFC data. This table includes the sum of the fixed and random coefficient estimates for each EFC. In general, there were slight but significant differences among coefficient estimates for the same scale but different EFCs, except for the slope for scale S4 and the intercept for the scale S5. These differences were even smaller for the S3 and S4 scales where the correction would be unnecessary. The greatest differences among EFCs and, therefore, the greatest improvements between marginal and conditional R 2 were found for scale S6.
Residuals for MSCcorr obtained after applying the bridging function, including both fixed and random effects, (Fig. 3b), can be compared with the residuals for MSCpred (Fig. 3a). The comparison showed that for large cover values estimated with scales S1 and S2, this correction would be important because it circumvents underestimation. For scale S5, it would be also important to apply the correction to the MSCpred values to keep the residuals within the limits of ± 5% and to avoid overestimation for large cover values.
However, for the scale S6, despite the improvement in results after correction, it can be seen that large errors were still found (MAE = 2.6 and RMSE = 3.4). In addition, the residuals were large mostly for the species less well represented in the EFC, while for dominant species whose occurrences were recorded for large numbers of plots, the residuals were less than 5% (Fig. 5).

Discussion
Forest monitoring has tended to focus on the tree component with little attention paid to understory components. However, ecosystem services and productivity for understory vegetation are probably comparable to those of the tree component (Nilsson and Wardle 2005). Additionally, because forest management practices that alter fire regimes and the composition of shrub vegetation may have long-term consequences for both conservation goals and commercial forest productivity (Nilsson and Wardle 2005), robust and updated information is necessary. In recent decades, the scopes of NFIs have broadened to include new variables (Tomppo et al. 2010). However, estimates produced by different countries are often not comparable because of differences in NFI definitions, plot configurations, measured variables, and measurement protocols. As a consequence, harmonizing estimates produced at national levels is essential for the production of sound EU forest information (Vidal et al. 2016b).
For this study, MSC was defined as a key forest indicator primarily because of its relevance, for being informative over large regions (Karl et al. 2017), but also because of the large number of countries with available information to estimate it. The proposed shrub indicator, MSC, can be estimated for species from NFI data for all except three of the countries that assess shrub cover. Additionally, two countries consider only a very few shrub species, so the use of the indicator would be restricted to those species. This indicator facilitates characterization of EFC understories and provides important ecological information. Additionally, it would be useful for forest management for estimating shrub biomass and carbon by recording additional variables such us average shrub height (Pasalodos-Tato et al. 2015).
Species monitored by individual countries are generally selected according to their inventory aims and the understory structure and composition of their forests . For example, Germany records forest plant species that cause forest management problems, while Spain considers both floristic and ecological aspects, as well as other factors of interest such as biomass, wild-fires or livestock browsing. Therefore, shrub composition harmonization, or even standardization, would require changing national-level field protocols. Because MSC estimates are difficult to compare unless the same species are assessed, a species list that is standardized across countries would be particularly useful.
As expected, the analyses of the scale effect for estimation of MSC revealed a relationship between the residuals and the widths of the cover scale intervals: The wider the scale intervals, the greater the residuals (Table 4). In addition, the different shapes of the empirical cover distribution curves affected the results: Downward concavity resulting from most plots having small species cover values led to the need for more accurate scales for the smaller ranges of cover (less than 50%) (Figs. 2, 3). Maximum values of species covers were generally less than 50%; therefore, use of S6 or S5 with fewer scale intervals in the range between 10 and 50% led to greater residuals (Table 4). In fact, MSC values for most of the species ranged between 10 and 25%. Due to the asymmetric distribution, the use of S5 and S6 with only one interval between 10 and 50, or 0 and 50%, overestimated the indicator values, whereas for the other scales they were generally underestimated (Fig. 4). This is also the reason that S3 and S4 produced smaller MEs and MAEs, and smaller residuals distribution, i.e., because they are the only scales having intervals in the range (5, 25%]. Nevertheless, MAE and RMSE indicate that mean errors for S1, S2, S3 and S4 were usually less than 3%, even when significant differences were found. This value can be considered acceptable when considering that visual error estimates due to observer bias vary greatly among species (Bergstedt et al. 2009;Kercher et al. 2003;Vittoz et al. 2010). Furthermore, a visual assessment with 1% interval widths is not realistic; for instance, Spanish NFI field teams tend to round to the nearest 5% when estimating cover (Figs. 1, 2, 3), an effect also noted by Wilson (2007).
Although interval midpoints for the cover scales have been commonly used as interval marks for transformation between scales (Canullo et al. 2011), other values representing each cover interval could be considered. Midpoints have not been demonstrated to provide the most accurate values. In fact, the results of this study revealed that maximum values of MSC and species covers are generally less than 25 and 50%, respectively, meaning that for the cover scales with only one interval up to 50% (e.g., S5 and S6), the interval mark should be less than the midpoint. For this reason, and with the objective of minimizing the scale effect, bridging functions were formulated as linear mixed models with the EFC as the grouping structure. Results showed that after application of bridges, MSCs were comparable for most of the NFI scales. S3 and S4 provided accurate MSC estimates even without applying the bridging functions. Moreover, the results from the validation process demonstrated that the coefficient estimates could be used for correction for independent datasets (Table 4, Fig. 3). Figure 5 shows that the MSC residuals following application of the bridging function decreased when the number of plots for which a species cover was recorded increased. However, when the number of plots was less than 20, residuals were greater than 5% for S5 and S6. Therefore, the results of Fig. 3 Mean species cover residuals for the six different analyzed scales, a before the bridging function correction and b after applying the bridging function. Mean error (ME), mean absolute error (MAE) and root mean square error (RMSE) are presented. a. Before bridging function. b. After bridging function ◂ Fig. 4 MSC residuals, i.e., observed minus predicted, means and confidence intervals at 95% for the six different analyzed scales by EFC this study could be applied to species with limited distribution areas and presented in few sample plots, but only when their covers were recorded using scales S1, S2, S3 or S4.
The study was conducted using only data from the Spanish NFI, which for shrub purposes uses plots with a 10-m radius. Therefore, it would be interesting to conduct further analyses using data from more countries and forest types. Particular emphasis in future analysis should focus on the effect of the size of NFI monitoring areas which vary considerably among European countries (Alberdi et al. 2010). This size effect was noted by Klimeš (2003) as one of the reasons for variation in visual cover estimates. Deviations between countries have also been noted in the Level I European network established in the frame of the International Co-operative Programme on the Assessment and Monitoring of Air Pollution Effects on Forests, ICP Forests (Van Dobben and De Vries 2010). This phenomenon is attributed to both use of different methods such as distinct cover scales and also the observer effect resulting from differences in taxonomic views and experience (Archaux et al. 2009).
Finally, to facilitate use of MSC for harmonized reporting, two standardization recommendations are proposed. First, standardized NFI shrub species lists with common criteria for the selection of species for these lists are needed. We propose using the percentage scale with 1% interval widths (S0), although scales S1, S, S3, S4 and S5 also exhibited small mean errors and residuals. Wilson (2007) also recommended against using cover classes because they contribute to loss of accuracy and are not much faster. Though visual assessments using a cover scale with intervals of 1% width rather inaccurate in relation to the scale, statistical analysis would be facilitated, and the errors at the interval margins would be reduced (Wilson 2007). Second, addition of a variable that facilitates estimation of shrub biomass is recommended; possibilities include average height by species

Conclusions
Mean species cover by EFCs can be used as an indicator for harmonized assessments. In the framework of European forest information, MSC could be estimated by most countries using data from their respective NFIs. The scale effect revealed that greatest residuals occurred in the smaller ranges of the cover scales (less than 50%) which meant that more and narrower scale intervals for the small classes yield smaller residuals. With the exception of scales S6 and S5, the differences between observed and predicted indicator values were less than 5%.
S3 and S4, the scales most similar to the Braun-Blanquet scale, could provide results for the indicator even without harmonization. The bridging function based on linear models improved the harmonization of MSC for the other scales, even for more unfavorable situations associated with scales having only two intervals (S6). The use of these bridging functions could be transferred to other countries although it would be recommendable to determine specific national bridging functions, especially for most unfavorable scales, i.e., S5 and S6.
Although this study is focused on NFIs that provide national assessments, the MSC indicator can also be used to harmonize estimates at local and regional scales. However, for species that are present in few sample plots, only cover scales S0, S1, S2, S3 and S4 can be considered.
Because models were obtained from an experimental approach using data from the SNFI, and because they could differ among countries and forest types, further analysis is recommended.
Acknowledgements The research was funded by the projects DIAB-OLO (European Union's Horizon 2020 research and innovation program under Grant Agreement No. 633464) and EG12-0073 (INIA). We wish to thank all the people who compiled the information on forest biodiversity in NFIs within the of COST Action E43 "Harmonisation of National Inventories in Europe: Techniques for Common Reporting" (Chairman Prof. Erkki Tomppo) within Working Group 3: "Harmonised indicators and estimation procedures for assessing components of biodiversity with NFI data" led by Gherardo Chirici. We also want to express our thanks to Roberto Vallejo and Vicente Sandoval of the Spanish Ministry of Agriculture, Food and Environment for kindly providing access to the full Spanish NFI datasets and to Paula Gil, Laura Hernández, Isabel Cañellas, Alfonso San Miguel and Miren del Río for their inestimable suggestions. Finally, the authors wish to Fig. 5 Relationship between the number of plots where a species were recorded and the residual for MSC after applying the bridging function, for the six different analyzed scales thank Roberto Canullo and Han van Dobben for their collaboration to improve this paper.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.