National Forest Inventories capture the multifunctionality of managed forests in Germany

Abstract

Background

Forests perform various important ecosystem functions that contribute to ecosystem services. In many parts of the world, forest management has shifted from a focus on timber production to multi-purpose forestry, combining timber production with the supply of other forest ecosystem services. However, it is unclear which forest types provide which ecosystem services and to what extent forests primarily managed for timber already supply multiple ecosystem services. Based on a comprehensive dataset collected across 150 forest plots in three regions differing in management intensity and species composition, we develop models to predict the potential supply of 13 ecosystem services. We use those models to assess the level of multifunctionality of managed forests at the national level using national forest inventory data.

Results

Looking at the potential supply of ecosystem services, we found trade-offs (e.g. between both bark beetle control or dung decomposition and both productivity or soil carbon stocks) as well as synergies (e.g. for temperature regulation, carbon storage and culturally interesting plants) across the 53 most dominant forest types in Germany. No single forest type provided all ecosystem services equally. Some ecosystem services showed comparable levels across forest types (e.g. decomposition or richness of saprotrophs), while others varied strongly, depending on forest structural attributes (e.g. phosphorous availability or cover of edible plants) or tree species composition (e.g. potential nitrification activity). Variability in potential supply of ecosystem services was only to a lesser extent driven by environmental conditions. However, the geographic variation in ecosystem function supply across Germany was closely linked with the distribution of main tree species.

Conclusions

Our results show that forest multifunctionality is limited to subsets of ecosystem services. The importance of tree species composition highlights that a lack of multifunctionality at the stand level can be compensated by managing forests at the landscape level, when stands of complementary forest types are combined. These results imply that multi-purpose forestry should be based on a variety of forest types requiring coordinated planning across larger spatial scales.

Background

Forests supply multiple regulating, material and non-material ecosystem services, i.e. nature’s contributions to people (IPBES 2019), including timber and other non-timber products, wild food, carbon sequestration, groundwater recharge, flood regulation, protection from soil erosion, recreational opportunities and habitat provision (Bauhus et al. 2010; Gamfeldt et al. 2013; Miura et al. 2015; Mori 2017; van der Plas et al. 2017; Storch et al. 2018). Both national (BMELV 2011) and international (FAO 2013) policy guidelines call for new strategies in forest management, i.e. multi-purpose forestry, to deliver as many of these ecosystem services as possible simultaneously. Yet, the major goal in managed forests is typically the production of timber and other wood products, which provide the main or only source of income for forestry. In addition, it is likely that not all ecosystem services can be maximized simultaneously, due to trade-offs between different services (van der Plas et al. 2016; Mouchet et al. 2017; van der Plas et al. 2017; Turkelboom et al. 2018). Nevertheless, recent studies found a large potential for forest multifunctionality at local (Felipe-Lucia et al. 2018) and continental scales (van der Plas et al. 2017). In order to unlock this potential, it is crucial to understand the degree to which managed forests already reflect multi-purpose forestry and how multifunctionality differs at the stand scale between different forest types.

Generally, multi-purpose forestry aims at the supply of additional ecosystem services (i.e. benefits humans obtain from ecosystems) together with timber production, such as non-timber products (e.g. berries, mushrooms, game) and recreational activities. The supply of ecosystem services can be measured either directly through demand for or use of an ecosystem service, or indirectly through ecosystem functions or processes which contribute to ecosystem services directly or indirectly (Garland et al. 2020). The level of multifunctionality of a forest can hence be assessed either as the potential of a forest to supply multiple ecosystem functions (stand scale or ecosystem-function multifunctionality) or as the actual use of the forest for multiple activities (landscape scale or ecosystem-service multifunctionality) (Manning et al. 2018). Here, we model the potential supply of ecosystem services independently of their actual use or demand (i.e. ecosystem-function multifunctionality) based on forest attributes such as stand density or tree species composition and other abiotic variables collected across a range of forest plots in Germany. Ecosystem services which are mainly driven by demand or (cascading) uses, such as tourism, long-term climate change mitigation or replacement of fossil fuels, were not included as data were not available at the relevant local scales. We use those models (also termed ‘production functions’ (Nelson et al. 2009; Tallis and Polasky 2011)) to predict the potential supply of forest ecosystems across Germany, using data from the National Forest Inventory (NFIs). NFIs provide standardized values for forest attributes and have long been used to assess the success of forest management strategies (Vidal et al. 2016a; Vidal et al. 2016b) and have been explored for assessments of ecosystem services or their underlying functions (Gamfeldt et al. 2013; Corona 2016; van der Plas et al. 2017; Storch et al. 2018). Based on the predicted potential supply of ecosystem services, we evaluate the multifunctionality of managed forests at the national level.

Methods

Calibration sites

We sampled 150 forest plots of 100 m × 100 m distributed across three regions in Germany (Schwäbische Alb in the South-West, Hainich-Dün in the Center and Schorfheide-Chorin in the North-East of Germany) and covering the dominating forest types in Central Europe up to 800 m elevation, except floodplain forests. The plots are part of the long-term research platform Biodiversity Exploratories (Fischer et al. 2010) and include a range of management intensities from unmanaged forests to conifer plantations.

German National Forest Inventory (predicted) sites

We used data from the most recent German National Forest Inventory (NFI) in 2012. The NFI was conducted on 27,121 plots arranged on a regular 4 km × 4 km basic grid. At each grid, an inventory cluster of four plots with a side length of 150 m is located, with the cluster coordinate indicating the location of the south-west plot. Each of the four plots is defined by a central sampling point for measurements on individual trees and circles of different radii for measurements of additional forests structures (Polley 2011). Only clusters with at least one plot located within a forested area were surveyed. For this study, we also excluded plots where access is hindered by difficult site conditions or restricted due to protected areas.

Potential supply of ecosystem services

In each of the calibration sites, we collected data on indicators (sensu Garland et al. 2020) of the potential supply of 15 ecosystem services (Supplementary Table 1), which cover the three main categories established by the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services IPBES (IPBES 2019) (regulating, material and non-material services). We included indicators for ‘climate regulation’ (carbon storage in trees, soil carbon stocks, local temperature regulation); ‘formation, protection and decontamination of soils and sediments’ (root decomposition, dung decomposition, potential nitrification activity, phosphorus availability, mycorrhiza and saprotrophic fungal richness); ‘regulation of detrimental organisms and biological processes’ (bark beetle control); ‘food and feed’ (edible fungal richness, cover of edible plants); ‘materials, companionship and labor’ (forest productivity, as a proxy for timber production); as well as ‘learning and inspiration’ (cover of culturally interesting plants, bird species richness). While the selected set of indicators represent the potential supply of multiple forest ecosystem services, our sampling procedures could not assess the actual use and demand of other (more sensitive) ecosystem services, such as hunting or recreational value.

Soil-related ecosystem services were sampled in a joint soil-sampling campaign that took place in May 2011. Within each of the 150 plots, 14 samples were taken of the upper 10 cm of the mineral soil along two 40 m long transects using cores with a diameter of 5 cm. The 14 soil samples were mixed into a composite sample before further analysis. Mycorrhizal and saprotrophic fungi richness (Buscot et al. 2018a, 2018b; Wubet et al. 2018; Schröter et al. 2019), richness of edible fungi (Allan et al. 2014), potential nitrification activity (Schloter and Stempfhuber 2018), phosphorus availability (Sorkau and Oelmann 2018), and soil carbon (Schöning et al. 2018a) were determined from those composite soil samples.

Richness of mycorrhiza and saprotrophic fungi

Saprotrophic fungi contribute to nutrient turnover by decomposing the organic material produced by plants (Baldrian and Valaskova 2008) while mycorrhizal fungi facilitate the plant nutrient uptake in return for photosynthesized carbon, which they channel into the soil, making a major contribution to long term C sequestration (Clemmensen et al. 2013; Pena 2016).

Higher species richness of soil fungal communities has been associated with higher functional diversity (Courty et al. 2010; Clemmensen et al. 2013) and therefore healthier forests, as fungi have different roles in decomposition (van der Wal et al. 2013), water and nutrient uptake (Courty et al. 2010; Pena et al. 2013; Pena and Polle 2014), enzyme production (Buee et al. 2007; Pritsch and Garbaye 2011), and soil health (Lehmann et al. 2019). We sampled the whole fungal community following the above-mentioned description of the joint soil-campaign campaign. DNA was separately extracted from soil and root samples; afterwards amplified using ITS-primer sets and prepared for 454 pyrosequencing (Goldmann et al. 2015; Schröter et al. 2015). To cope with different sequencing approaches bioinformatically, raw root and soil sequences were initially trimmed separately using MOTHUR (Schloss et al. 2009): ambiguous bases, homo-polymers and primer differences of more than eight bases were removed; all primer and barcode sequences were discarded; and sequence reads with a quality score lower than 20 and read length less than 360 bp were removed. Afterwards, ITSx was used to identify, cut and align the ITS2 region in both sequence sets (Bengtsson-Palme et al. 2013). Moreover, sequences identified as plant-borne were removed. This was important, since plant contaminations particularly occurred in sequencing of root fungi. Then, we performed a chimera check of both datasets using the uchime algorithm (Edgar et al. 2011) implemented in MOTHUR. Potential chimeric sequences were discarded subsequently. In order to gain shared operational taxonomic units (OTU) from root and soil sequences, the two datasets were combined before the clustering. VSearch (Rognes et al. 2016) was executed to obtain OTUs using 97% sequence similarity. An additional chimera check was carried out thereafter (again using the uchime algorithm implemented in MOTHUR). The taxonomical assignment was done using MOTHUR against the UNITE fungal database (version 7.2) (Kõljalg et al. 2013). In addition, the database FunGuild (version 1.0) (Nguyen et al. 2016) was used for functional assignment of fungal OTU. The richness of saprotrophic fungi OTUs and mycorrhizal fungi OTUs were used to estimate the potential supply of ecosystem services.

Richness of edible fungi

Mushroom collection or observation is common in forests, and in addition to providing a type of wild food, it is an important recreational activity (Millennium Ecosystem Assessment 2005; Haines-Young and Potschin 2013). We estimated the potential of our forests to harbor edible fungi by analyzing fungal species pools in forest soils following the abovementioned description of the joint soil-sampling campaign. Edible fungi were identified following the criteria of the German Mycological Society, excluding those species with inconsistent edible value (Deutsche Gesellschaft für Mykologie e.V. 2015). We used species richness of edible fungi as a proxy of potential edible fungi observation (See complete list in Felipe-Lucia et al. (2018) Supplementary Table 9).

Potential nitrification activity

Nitrogen (N) can be a limiting element for plant growing and therefore limit the functioning of the ecosystem (Vitousek and Howarth 1991; LeBauer and Treseder 2008). We investigated the nitrification process in forest soils in terms of potential nitrification activity. Soil samples were collected during the joint soil-sampling campaign as described above. Potential nitrification measures were derived from the abundance of nitrifying bacteria following (Hoffmann et al. 2007) and used as a proxy for potential nitrification activity (Allan et al. 2015; Soliveres et al. 2016). Two of the readings were excluded from the analysis because the standard deviation was three times larger than the mean and hence considered unexplained outliers.

Available phosphorus

Due to continuously high atmospheric N deposition, phosphorus (P) becomes increasingly important as a limiting element for plant growth and therefore might limit the functioning of the ecosystem (Holland et al. 2005; Vitousek et al. 2010; Lang et al. 2017; Clausing et al. 2020). We investigated the availability of P in forest soils by collecting soil samples as described above. Total P was extracted with 0.5 mol∙L^− 1 NaHCO₃ (pH = 8.5) following the Olsen methodology (Olsen 1954; Alt et al. 2011; Sorkau et al. 2018) and measured using Inductively Coupled Plasma/Optical Emission and Spectrometry (ICP-OES, PerkinElmer Optima 5300 DV, S10 auto sampler). P concentrations in the extraction solution was used as a proxy of P availability for plants (Felipe-Lucia et al. 2014).

Carbon stocks in the soil

Forest soils are important carbon pools. In order to estimate the amount of carbon stored in forest topsoils, we followed the abovementioned joint soil-sampling campaign. An aliquot of < 2 mm sieved soil was homogenized with a ball mill (RETSCH MM200, Retsch, Haan, Germany) and used to determine total C concentrations by dry combustion in an elemental analyzer (VarioMax, Hanau, Germany). Inorganic carbon was determined after combustion of organic carbon at 450 °C for 16 h using the same elemental analyzer and total organic carbon was afterwards determined from the difference between total and inorganic carbon (Schöning et al. 2018a). Organic carbon stocks were determined by multiplying organic carbon concentrations with the total soil mass (0–10 cm) per unit area (< 2 mm) per m² in each plot.

Root decomposition

Fine root decomposition plays an important role in element cycling in forest ecosystems (Hobbie 1992). We measured decomposition of fine roots (< 2 mm) within the upper 10 cm of the mineral soil. Three polyester litterbags per plot, with a mesh size of 100 μm, were filled with fine roots collected from 2-year-old European beech (Fagus sylvatica L.). These were buried in each of the plots in October 2011 and were then harvested after 12 months in October 2012. We used the percentage of root litter mass loss as a proxy of decomposition (Solly et al. 2011; Solly et al. 2014).

Dung decomposition

Dung beetle communities contribute to the rapid decomposition of fecal deposits from both wild mammals and domestic livestock, representing a key ecosystem service. We installed five dung piles (cow, sheep, horse, wild boar, red deer) on each plot and collected the remaining dung after 48 h. The average percentage of dung dry mass removed (mostly by tunneling dung beetles) was used as an indicator of dung removal rates (Frank et al. 2017; Frank and Blüthgen 2018).

Carbon storage in living trees

Trees are important carbon sinks as they store carbon in their tissues via photosynthesis. In order to assess the amount of carbon stored in trees, we estimated the living tree volume on each plot as assessed in the second forest inventory (2015–2016) (Kahl and Bauhus 2014; Schall and Ammer 2017). The living tree volume on each plot which was converted to dry biomass using standard conversion factors of 0.46 for conifers (average of spruce (0.43) and pine (0.49)) and 0.67 for broadleaved trees (average of oak (0.66) and beech (0.68)), according to Lohmann (2011). The carbon stored in a tree is approximately 50% of its dry biomass (The Intergovernmental Panel on Climate Change (IPCC) 2003).

Bark beetle control

Natural bark beetle control is an important forest ecosystem service that can have an effect on other services like production of quality timber and aesthetic value (Jactel et al. 2009; Bengtsson 2015). We assessed the abundance of potential pest species among the bark beetles (i.e. ambrosia beetles), together with their antagonists. Ambrosia beetles cause substantial damage worldwide are thus considered as important pest species (Grégoire et al. 2015). In Europe they impact broad-leaved (e.g. Xyleborus dispar, Xylosandrus germanus, Trypodendron domesticum) as well as conifer trees (e.g. Trypodendron lineatum). Bark beetle control was estimated from the ratio of bark beetle antagonists to bark beetles (i.e. tribe Xyleborini) based on collections of bark beetles and their antagonists in pheromone traps (Weisser and Gossner 2017). The collected specimen included species which attack broad-leaf trees as well as species which attack conifer trees (Gossner et al. 2019). Lineatin lures and ethanol were used as attractants for bark beetles and their antagonists. Traps were emptied every second day during the main activity period of bark beetles in 2010 (Gossner et al. 2019) and in weekly to monthly intervals afterwards. We standardized the data based on the method proposed by Grégoire et al. (2001). We only considered bark beetles that are attracted by lineatin and/or ethanol (i.e. species within the tribe Xyleborini) and predators and parasitoids that are mentioned as antagonists of Xyleborini in the literature (Kenis et al. 2004). Predator-prey ratios have been frequently used as measure of pest control potential in different systems (van der Werf et al. 1994; Klein et al. 2002; Bianchi et al. 2013). Bark beetles vs. predators and parasitoids ratio was used as a proxy of bark beetle control.

Temperature regulation

Forests buffer extreme temperatures due to their dense canopy cover (Frey et al. 2016). We collected data on air temperature 2 m above ground from climatic stations installed in each of the 150 plots. Records were quality controlled with respect to temporal dynamics and by using data from neighboring stations. This results in data for 143 plots suitable for the analysis (Hänsel and Nauss 2019). Temperature regulation was defined as the inverse of the diurnal temperature ranges (DTR) and calculated as the difference between daily maximum and minimum temperature values (Scheitlin and Dixon 2010). The inverse value of the average DTR per plot (1/DTR) was used as the proxy in order to facilitate the interpretation of the results, so higher values of the proxy mean higher temperature regulation of the forest plot (Felipe-Lucia et al. 2014). We used only 2011 as most other indicators of ecosystem services were assessed in 2011.

Forest productivity

Timber is one of the main products extracted from forests. Since accurate harvesting records were not available for our calibration sites, we used the increment of the stand in wood volume per hectare and year as a proxy for forest productivity. We quantified productivity for each of the 150 plots as mean annual increment (MAI) across rotation (i.e. culmination of MAI) for even-aged forests and as periodic annual increment (PAI) between two forest inventories for uneven-aged and unmanaged forests (Schall and Ammer 2018a, 2018b). MAI was estimated based on site class or site maps of forest administrations. Culmination of MAI is estimated on 70 to 100 years for Norway spruce (Picea abies), 70 to 90 years for Scots pine (Pinus sylvestris), 120 to 140 years for oak (Quercus spp.), and 140 to 160 years for European beech (Fagus sylvatica). PAI was estimated as the difference between the increment measured during the first forest inventory (2008–2011) and the second forest inventory (2015–2016) of our plots divided by the time span in years. All values are given as volume above bark (> 7 cm in diameter) in m³ per ha and year.

Edible vascular plants

A common recreational use of forests is the collection of fruits, nuts, berries and other plant parts for cooking (Gamfeldt et al. 2013). We estimated the potential supply of edible vascular plants based on two vegetation surveys (in spring and summer to represent both flowering aspects) conducted in a 20 m × 20 m subplot in each of the plots in 2011 (Schäfer et al. 2017). Cover of each plant species was estimated across four different vegetation layers (herbs, shrubs < 5 m, trees 5–10 m, trees > 10 m), and the cumulative cover per vegetation layer was used to calculate overall cover for each species. Wild edible plant species known to be collected were identified by botanists from the Botanical Society of Bern (Bernische Botanische Gesellschaft) with knowledge on people preferences. Total cover of these species was used as a proxy of potential wild edible plants gathering (See complete list in Felipe-Lucia et al. (2018) Supplementary Table 10).

Culturally interesting vascular plants

Plants blooming in early spring in forests are highly appreciated for their aesthetic value (Bhattacharya et al. 2005). Other plant species are of special interest for interested non-experts and botanists, such as the forest specialists Helleborus spp., Asarum europaeum, Gallium odoratum or Gagea lutea, Hepatica noblis and Anemone nemorosa because of their unique flowering times or their medicinal use (Schmidt et al. 2011; Lauber et al. 2012). We estimated the potential aesthetic and educational value of forest plots following the abovementioned methods description for edible vascular plants. Plant species of special interest for the general public or for botanists were identified by botanists from the Botanical Society of Bern (Bernische Botanische Gesellschaft) with knowledge on people preferences. Total cover of these species was used as a proxy of potential aesthetic value (complete list in Felipe-Lucia et al. (2018) Supplementary Table 11).

Bird species richness

Forests provide good opportunities for birdwatching. In order to estimate the bird-watching potential of forest plots, we performed five bird surveys during breeding times (March to June 2011; Renner and Tschapka 2017) counting the number of individuals seen or heard during 5 min at the center of each plot (Renner et al. 2014). The number of bird species was used as a proxy for bird watching potential.

Forest structural and compositional parameters and abiotic parameters

We selected 18 forest parameters shown to affect ecosystem functions and services that were available for both the calibration sites (Biodiversity Exploratories sites) and the predicted sites (German NFI plots). The parameters are related to stand structure (stand density, aggregated values of diameter at breast height, diversity of diameter classes, basal area of trees, relative crown cover) (Schall and Ammer 2017), tree species composition (number of tree species, relative cover of conifers, evenness of tree species) (Schall and Ammer 2017), regeneration (number and evenness of woody species and young trees, relative cover of young Fagus sylvatica (L.) in the understorey) (Grassein and Fischer 2013), and deadwood (volume of coarse deadwood, diversity of deadwood types) (Kahl and Bauhus 2018).

All forest parameters were calculated on the calibration sites based on a full inventory of all trees with a diameter at breast height (dbh) larger than 7 cm within the 1-ha plots. In contrast, trees were assessed using angle-count sampling or within circles of different size depending on the height of the tree on the predicted sites (Schmitz et al. 2008). Hence, not all trees within 1 ha were assessed and measured for the predicted sites and parameters had to be extrapolated to 1 ha from the individual tree counts. From each individual tree count, the density per hectare of similar trees was calculated from the angle-counts. Extrapolation based on an average of four angle-counts gives very good estimations of parameters of 1 ha, except for parameters which are based on the diversity of trees in size and species identity (Supplementary Fig. 1). For those parameters, the sampling effect (higher diversity with more sampled individuals) leads to much higher values from the full inventory compared to the angle-count assessment. Hence, we used the average value of four simulated angle-counts (instead of the full inventory) on the calibration sites for forest parameters related to diversity.

The stand density of all trees and of the main tree species is measured as the number of tree individuals per hectare which have a diameter at breast height (dbh) (at 1.3 m) larger than 7 cm. The following aggregated values of the dbh were calculated: arithmetic mean, standard deviation, quadratic mean and arithmetic mean among the 50 largest trees. The standard deviation and coefficient of variation were estimated based on the virtual angle-count method.

Staudhammer and LeMay (2001) recommend the use of a tree size diversity index which is calculated as the arithmetic mean of the Shannon diversity (H′) based on height classes (H′_h) and the Shannon diversity based on diameter classes (H′_d). However, as information on tree height was not available for 22 plots (mostly thickets with high tree densities), we only used the H′_d based on classes of diameter at breast height (at 1.3 m) in steps of 4 cm (beginning with 7 cm). The basal area (m²) per hectare of all trees within a size class was used as the abundance for the Shannon diversity. The diversity of diameter classes was estimated based on the virtual angle-count method. The number of tree species was counted based on trees from the virtual angle-count method. The overall crown projection area was calculated by multiplying each tree’s crown projection area with its density per hectare. The relative crown projection area was calculated relative to 1 ha. The relative cover of conifers is based on the covered ground by all trees, which is calculated for a normalized reference area (‘ideeller Flächenanteil’). Coniferous species in the Biodiversity Exploratories mostly include spruce and pine. In the German National Forest Inventory, the following coniferous species are listed: Norway spruce (Picea abies), Scots pine (Pinus sylvestris), Sitka spruce (Picea sitchensis), Douglas fir (Pseudotsuga menziesii), Larches (Larix spec.), Yew (Taxus baccata), other spruces (Picea spec.), other pines (Pinus spec.), and other firs (Abies spec.). The evenness of trees was calculated as Pielou’s evenness. This measure is the Shannon diversity (H′) divided by the logarithm of species richness. The Shannon diversity was calculated based on tree density as a substitute for abundance. Evenness of trees was calculated based on the full inventory as extrapolations based on virtual angle-count measures were not biased (Supplementary Fig. 1). The number of woody species in the shrub layer was assessed within a square of 400 m² for the Biodiversity Exploratories. In the German NFI, the number of woody species in the shrub layer was based on the species assessed in a 10-m circle for ground vegetation assessment, which includes trees up to 4-m height. The evenness of woody species was calculated as Pielou’s evenness with Shannon diversity based on the percentage cover as a substitute for abundance. The cover of young European beech trees within the shrub layer was measured relative to the plot size, i.e. not relative to other species in the shrub layer. In the Biodiversity Exploratories, the shrub layer was assessed within a square of 400 m². In the German National Forest Inventory, the shrub layer is assessed within a 10-m circle. Only individuals of Fagus sylvatica which are below 4 m in height were considered. Total volume of coarse deadwood was estimated based on lying and standing deadwood items with a diameter of at least 25 cm. Smaller deadwood items and tree stumps were not considered. The volume was calculated based on the length or height of each deadwood item and its diameter. The German National Forest Inventory uses the assumption that the deadwood item resembles a cylinder, thereby averaging between thinner or thicker parts of the deadwood item (Schmitz et al. 2008):

$$ \mathrm{Volume}\ \left({\mathrm{m}}^3\right)=\frac{\uppi \times {\left[\mathrm{diameter}\ \left(\mathrm{cm}\right)\right]}^2\times \mathrm{length}\ \left(\mathrm{dm}\right)}{4\times {100}^2\times 10} $$

As all deadwood items on the plot were assessed in the Biodiversity Exploratories, the sum of all deadwood volumes indicates the deadwood volume per hectare. As deadwood items are only assessed within a circle of 5-m radius in the German National Forest Inventory, the volume per hectare must be calculated with a correction factor:

$$ \mathrm{Volume}\ \left(\frac{{\mathrm{m}}^3}{\mathrm{ha}}\right)=\frac{\mathrm{volume}\ \left({\mathrm{m}}^3\right)\times 1000}{F_{500}^{korr}}\ \mathrm{with}\ {F}_{500}^{korr}=78.54\ {\mathrm{m}}^2 $$

Diversity of deadwood was estimated as the Simpson diversity based on types of deadwood (standing, downed tree, downed stem, other deadwood), the functional tree group (coniferous or deciduous) and decay stage (not decayed, beginning decay, proceeding decay and heavily decomposed). In the German National Forest Inventory, 40 unique combinations of categories were found, the Biodiversity Exploratories had 33 unique combinations of categories. Each combination is treated as a species and its overall volume is treated as abundance to calculate the Simpson diversity.

Abiotic parameters

Abiotic variables included elevation above sea level, aspect and slope at the center of both the calibration and the predicted sites (Nieschulze et al. 2018a; Nieschulze et al. 2018b). In addition, topsoil properties (sand, silt and clay content) as well as soil depth were used (Schöning et al. 2018b). For the Biodiversity Exploratories, elevation above sea level (m) was derived from a digital terrain model and taken at the center of the 1-ha plot. For the German National Forest Inventory, the inventory plot coordinates were projected onto a digital terrain model of Germany (Geodatenzentrum, http://www.bkg.bund.de, downloaded 14th August 2017) with a grid size of 200 m from which the elevation above sea level was extracted for the center of the plot. The aspect is given as the circular average over the plot area with 360° equaling 0° as due north. The slope is given in percent. For the Biodiversity Exploratories, aspect and slope were derived from a digital terrain model. Second-order finite differences were used to derive slope and aspect, all values are means calculated over the entire plot area. For the German National Forest Inventory, the inventory plot coordinates were projected onto a digital terrain model of Germany (Geodatenzentrum, http://www.bkg.bund.de, downloaded 14th August 2017) with a grid size of 200 m from which the aspect and slope were extracted for the center of the plot. The topsoil properties are described by the sand (2–0.063 mm), silt (0.063–0.002 mm) and clay (< 0.002 mm) content in percent. For the Biodiversity Exploratories, topsoil properties were measured in the upper 10 cm of the mineral soil taken during a standardized joint soil-sampling campaign in 2011. For the German National Forest Inventory, the inventory plot coordinates were projected onto the map of topsoil physical properties for Europe based on LUCAS topsoil data¹⁰⁸ with a grid size of 500 m. The soil depth is measured in cm and describes the rooting capacity of the soils. For the Biodiversity Exploratories, we determined soil depth by sampling a soil core in the center of all plots in 2008. We used a motor-driven soil-column cylinder with a diameter of 8.3 cm for the soil sampling (Eijkelkamp, Giesbeek, The Netherlands). The combined thickness of all topsoil and subsoil horizons was used as a proxy of soil depth. For the German National Forest Inventory, the inventory plot coordinates were projected onto a map of soil depth (Bug 2015) for Germany which is derived from profile data on the land use stratified soil map of Germany at scale 1:1,000,000. The lower limit of a soil is bedrock, or a groundwater influenced horizon.

Selection of predictors and model development

All predictors showed normal distribution of residuals. To reduce multi-collinearity among the initial set of 18 forest structural and compositional parameters and 7 abiotic parameters, we applied a variance inflation analysis (Zuur et al. 2009) using the package ‘fsmb’ (Nakazawa 2017) in R v.3.3.0 (R Core Team 2016), removing the predictor with the highest variance inflation factor until all variance inflation factors were below a threshold of 3 (Zuur et al. 2009). After this procedure, 10 forest structural and compositional parameters and 4 abiotic parameters remained as predictors. Complete information on the final set of predictors was available for 124 of the 150 plots.

Our ecosystem service models were based on Generalized Linear Models. In order to ensure normal distribution of residuals, some indicators of ecosystem services were transformed prior to model fitting: phosphorus availability, potential nitrification activity, bark beetle control, and cover of culturally interesting plants were square-root-transformed; cover of edible plants was log-transformed. We applied model simplification using the stepAIC function to find the most parsimonious model for each ecosystem service indicator. Models did not include any interactions between predictors as we wanted to derive purely additive models. We calculated a set of metrics for model quality: R², RMSE (root mean squared error), residual standard error of the whole model, and standard errors of estimated coefficients. Two of our models (edible fungi richness and bird species richness) showed low R² values (less than 20% of explained variance) and were discarded from further analyses as key drivers of those ecosystem functions were missing from our set of predictors. The other 13 ecosystem service models had at least 30% of their variance explained and were retained. The additional model quality metrics confirmed this selection (Supplementary Table 2).

While model quality metrics, like residual errors, are used to assess model quality for explanatory models (i.e. models to understand mechanisms), estimation of model quality for predictions requires a comparison of observed and predicted values. We applied this by splitting the Biodiversity Exploratories data into training data and test data. The training data comprised a random subset of 70% of the plots, equally distributed across the three regions. The test dataset comprised the remaining plots. For each ecosystem service indicator, models with the initial set of 25 predictors were calculated using the stepwise selection procedure. The resulting coefficients were used to predict the potential supply of ecosystem services for the test plots. Average R² values of 10-fold cross-validations from the function ‘train’ in the R package ‘caret’ v6.0 were then used to assess the model’s predictive accuracy. The whole process was repeated 999 times with different sets of plots in the training and test dataset. For all ecosystem service indicators, the predictive accuracy of the models was equal or higher than 30% and considered adequate (Supplementary Table 2). However, the predicted values for edible fungi and bird species richness were consistently higher than the observed values (Supplementary Fig. 2), so we excluded them from predictions at the national level. We provide the coefficients of these ‘low quality’ models as a potential starting point for the development of better models but advise against their direct use.

Definition of forest types for predictions

Each plot of the NFI should be considered as one of many sampling points, i.e. an individual inventory point is not representative of the forest stand or forest type around it. Hence, forest structural and compositional parameters have to be aggregated across several inventory points in order to accurately describe the characteristics of a forest (Riedel et al. 2017). The set of inventory points over which parameters are aggregated are usually defined as forest types which can be based on data derived from the inventory or based on additional information. We defined forest types based on the large-scale environmental characteristics of pre-defined geographic areas (also termed growth areas) (vTI Agriculture and Forestry Research 2012), and refined those based on tree species composition. The geographic areas were defined based on their geomorphological characteristics (i.e. bedrock type and topography), climate and landscape history (vTI Agriculture and Forestry Research 2012). Of the 27,121 initial inventory plots, only 25,246 inventory plots could be associated with maps of abiotic information (i.e. soil characteristics and topography), hence only those were used to define forest types.

The German NFI covers both common and rare forest habitat types, but only three of the rare forest habitat types are considered to be adequately represented: Luzulo-Fagetum beech forests, Asperulo-Fagetum beech forests and Galio-Carpinetum oak-hornbeam forests (Riedel et al. 2017). Hence, all forest inventory plots located within other rare forest habitat types were excluded (23,304 plots remaining). We only included geographic areas in which at least 40 inventory plots had been assessed (23,190 plots within 52 areas). Tree species composition was defined based on the identity of the main tree species (Fagus sylvatica, Picea abies, Pinus sylvestris, Quercus spp. and others) and the functional tree groups among the admixed tree species (pure stand, only deciduous species, only coniferous species, deciduous and coniferous species). Of the resulting 2652 possible combinations of geographic areas and tree species composition, 1668 combinations were found among the forest inventory plots and 148 combinations were represented by at least 40 inventory plots (14,645 plots total). Of those 148 combinations, 115 included a tree species combination which was also found within the calibration sites (11,762 plots in total). Among those 115 combinations, 53 were represented by at least 40 inventory plots with complete information on any predictor variable (Supplementary Table 3). These 53 forest types cover 7003 inventory plots with complete information on forest and abiotic parameters (Supplementary Fig. 3). Those plots were used to calculate average values and standard errors of all forest and abiotic parameters per forest type.

Prediction of potential ecosystem service supply in Germany

The average value and confidence intervals of each structural attribute were calculated from individual NFI plots and the upper and lower 95% confidence intervals of all forest and abiotic parameters calculated for the 53 forest types were compared to the observed ranges within the calibration sites to check for outliers (Supplementary Fig. 4). The ranges of predictors across all forest types fell within the observed range of the calibration sites except for soil depth for seven forest types. Hence, we are confident that our extrapolations of the potential supply of ecosystem services will not be confounded by exceeding the observed range of predictors. In addition, both sets of plots were compared with a principal component analysis (PCA) to check for distinct differences in multi-dimensional space (Supplementary Fig. 5). The two datasets show considerable overlap in multi-dimensional space, indicating that the Biodiversity Exploratories cover a large proportion of forest structural combinations found among the 53 most dominant forest types in Germany. However, the inherent differences in field assessments leads to a larger range of understorey diversity and evenness values within the Biodiversity Exploratories data compared to the forest inventory data (PC2 axis). On the other hand, the forest inventory plots show a larger range in tree species richness and tree size variability (PC3 axis). Overall, our results should be considered conservative because of the strict selection procedure regarding forest types and number of inventory plots included to ensure reliability of the predicted potential ecosystem service supply. While a broader definition of forest types will likely increase spatial coverage, it will also increase the variability in predictors within forest types and hence the uncertainty of predicted potential ecosystem service supply.

Based on our ecosystem service models, we predicted the potential supply of ecosystem services for the 53 forest types. We used the R function ‘predict’ and back transformed those ecosystem service indicators which were normalized prior to modelling. For each forest type, we calculated the potential supply of each ecosystem service based on the average values of forest and abiotic parameters across the inventory plots within each forest type, resulting in 53 average levels of potential ecosystem service supply. To identify groups of ecosystem services delivered by the same forest type, i.e. synergies within forest types, we conducted a principal component analysis (PCA) across the 53 forest types. To map the levels of potential ecosystem service supply, the average values were assigned to the 7426 individual German inventory plots which are located within the 53 forest types (including plots with incomplete information on structural and compositional parameters). To compare levels of potential ecosystem service supply for each forest type in relation to the other forest types, we calculated the relative level of potential ecosystem service supply based on the maximum predicted level across forest types. For each forest type, we then counted the number of ecosystem service indicators which had a predicted relative level of 0.5, 0.75 and 0.9 of the maximum predicted level. To assess the influence of variability in forest and abiotic parameters on the predicted level of potential ecosystem service supply per forest types, we predicted the potential ecosystem service supply for individual inventory plots and calculated average and lower as well as upper confidence intervals for each forest type. To assess which forest types showed the highest level of multifunctionality, we counted the number of ecosystem service indicators with at least average levels of potential supply by a forest type. We used the locations of the forest inventory points to assess whether the current state of forests at the landscape-scale already fulfills the expectations of multi-purpose forestry for Germany.

Results

The ecosystem service models differed markedly between the ecosystem service indicators, both in terms of the predictors selected as important, and in terms of their direction and strength of effects (Fig. 1 & Supplementary Table 2). For example, relative crown projection area had positive effects on carbon storage in trees and temperature regulation, but negative effects on nitrification potential and cover of edible plants. Proportion of silt in the soil, soil depth, relative conifer cover and variability in tree diameter had a significant effect on the majority of ecosystem service indicators. Relative cover of beech and evenness of woody species in the understorey as well as diversity of deadwood types had significant effects on one or two ecosystem service indicators.

Fig. 1

Coefficients of ecosystem service models based on forest structural and compositional parameters and abiotic parameters. Points represent mean estimates, lines represent the 95% confidence intervals based on t-distributions (all variables were normally distributed). Non-significant (i.e. p > 0.05) estimates are shown in grey, parameters without points were eliminated from the respective models through stepwise model simplification. Intercepts are shown separately from parameter coefficients for clarity; note the different ranges of the x-axes. Note that models for edible fungi and bird richness explain very little of the variability in the data (Supplementary Table 2) and were not used to predict potential ecosystem services supply

Trade-offs and synergies

The first PCA axis for potential ecosystem service supply across 53 forest types in Germany explained 46.6% of the variability and revealed synergies between dung decomposition and phosphorus availability, and to a lesser extent with bark beetle control (Fig. 2). Those three ecosystem service indicators showed trade-offs with mycorrhiza richness, soil carbon stocks and to a lesser extent with forest productivity and edible plants. The second axis (which explained 25.2% of the variability) revealed strong synergies between temperature regulation, carbon storage in trees, culturally interesting plants and nitrification potential. Forest types were clearly clustered in the ordination depending on tree species composition, with Norway spruce (Picea abies (L.) H. Karst) and Scots pine (Pinus sylvestris L.) separating forest types along the first PC (principal component) axis, while broadleaf forests and conifer forests separated along the second PC axis (Fig. 2a). Forests dominated by Norway spruce were positively associated with root decomposition, edible plants and saprotrophic fungal richness (cf. Awad et al. 2019), while forests dominated by European beech (Fagus sylvatica L.) were positively associated with local temperature regulation, culturally interesting plants (e.g. Anemone nemorosa L.), carbon storage in trees and potential nitrification activity. The association between forest types dominated by Norway spruce and edible plant cover likely indicates high (historic) management intensity, as many edible plants prefer the semi-open areas originating from management disturbances (e.g. Rubus spec.) or the acidic soils usually found in conifer forests (e.g. Vaccinium myrtillus L.). Forest types did not show clear clustering based on the geographic area (Fig. 2b), indicating that large-scale environmental conditions (geomorphology, climate and landscape history) were less important for differences in levels of potential ecosystem service supply than tree species composition. Exceptions were the forest types from the northern and eastern lowland (areas 13 & 22), which formed a distinct cluster. Those areas are characterized by sandy soils, which only allow for low soil carbon stocks but achieve high phosphorous availability (Grüneberg et al. 2013). This particular soil texture also makes those areas suitable for growing Scots pine, explaining the strong association between pine-dominated forest types and phosphorus availability.

Fig. 2

Trade-offs and synergies between ecosystem service indicators. Principal components analysis based on 13 ecosystem service indicators associated with tree species composition a and geographic area b and predicted for the 53 most dominant forest types in Germany. Spruce = Picea abies ((L.) H. Karst), Pine = Pinus sylvestris (L.), Beech = Fagus sylvatica (L.), Oak = Quercus spp., D= deciduous tree species, C = coniferous tree species

Level of potential ecosystem service supply

Forest types did not only differ in terms of the main ecosystem services they supported, they also differed greatly in the level at which each service could be supplied (Fig. 3). For potential nitrification activity, phosphorus availability, bark beetle control and cover of edible plants, most forest types supplied only up to half of the level than the ‘best’ forest types did. For those services, a particular combination of tree species composition and large-scale environmental conditions (reflected by the geographic area) seems to be required to achieve high levels of potential supply. For some ecosystem service indicators (e.g. potential nitrification activity or edible plants), different forest types from the same geographic area showed strong differences in the level at which each service was provided (Supplementary Fig. 3), indicating the importance of tree species composition. For other ecosystem service indicators (e.g. decomposition, richness of saprotrophs, carbon storage in trees and temperature regulation), differences in the predicted levels between forest types were less pronounced (Fig. 3 & Supplementary Fig. 6).

Fig. 3

Relative level of potential ecosystem service supply for each of the 53 most dominant forest types in Germany. Values are given relative to the maximum predicted level across forest types within each ecosystem service indicator. Colors indicate tree species composition and points are jittered horizontally to show the different geographic areas. Spruce = Picea abies ((L.) H. Karst), Pine = Pinus sylvestris (L.), Beech = Fagus sylvatica (L.), Oak = Quercus spp., D = deciduous tree species, C = coniferous tree species

No single forest type was superior to others in providing above-average levels of potential supply across all ecosystem service indicators (Fig. 4). Nevertheless, as each forest type achieved high levels for a subset of ecosystem service indicators, some forest types are better suited than others to provide a specific set of ecosystem services at the stand scale.

Fig. 4

Number of ecosystem service indicators which are supplied at or above a given threshold by each of the 53 dominant forest types in Germany. In total, 13 ecosystem services were modeled. Thresholds of relative ecosystem service supply were calculated relative to the maximum predicted level across all forest types within a service indicator. Colors indicate tree species composition and points are jittered horizontally to show the different geographic areas. Spruce = Picea abies ((L.) H. Karst, Pine = Pinus sylvestris (L.), Beech = Fagus sylvatica (L.), Oak = Quercus spp., D = deciduous tree species, C = coniferous tree species

Multifunctional forest areas

As the different forest types in this study represent a range of geographic areas, the differences among forest types translated into clear regional differences in predicted levels of potential ecosystem service supply (Fig. 5). Forest productivity levels (based on biomass increment) were highest in Northern Germany and this region also showed highest levels of regulating services (root decomposition, potential nitrification activity, mycorrhiza richness and carbon in living trees). Central Germany showed high levels of regulating services as well, but for a slightly different set of indicators (dung decomposition, phosphorus availability, saprotrophic fungal richness and soil carbon stocks). Southern Germany showed highest levels of yet another set of regulating services (phosphorus availability, bark beetle control and temperature regulation). Among those three regions, Southern Germany showed a more heterogeneous pattern of ecosystem services supply compared to Central or Northern Germany with all ecosystem service indicators showing average to high levels of supply at least in some parts of Southern Germany (Fig. 5b). These patterns correspond to the geographic distribution of the different forest types, with Southern Germany showing a larger variety of forest types compared to the Center and North (Fig. 5a). Hence, multi-purpose forest landscapes of a few hundred square-kilometers already exist in parts of Germany, and could be achieved in other parts as well by increasing the mix of tree species composition (under the constraints of geographic area and environmental conditions).

Fig. 5

Average predicted levels of ecosystem service indicators for 7426 forest inventory plots among the 53 dominant forest types in Germany. Upper maps a show distribution of tree species composition and geographic areas, which were used to define the forest types. Other maps b show the relative levels of each ecosystem service indicator; grey points indicate forested areas for which predictions were not possible (see Supplementary Fig. 3). Spruce = Picea abies ((L.) H. Karst, Pine = Pinus sylvestris (L.), Beech = Fagus sylvatica (L.), Oak = Quercus spp., D = deciduous tree species, C = coniferous tree species

Discussion

No single forest type supplied all studied ecosystem services at high levels and the specific combination of which services could be supplied depended on the forest type. To achieve high levels of multifunctionality at the regional or landscape-level scale (i.e. a spatial scale that includes numerous forest stands), multiple-purpose forestry needs to consider a mix of forest types, particularly of different tree species compositions in order to provide a wide range of ecosystem services.

The fact that each ecosystem service indicator is positively affected by a unique set of forest attributes and abiotic variables further indicates that multifunctionality within forest types can only be achieved for subsets of ecosystem services. Although tree species composition and site conditions are not independent, tree species composition was more important in determining the potential supply of ecosystem services than large-scale environmental conditions. A change in forest type or tree species composition will likely not result in a loss of supply for all ecosystem services. However, some ecosystem service indicators showed substantial variation within forest types, highlighting that variability in forest attributes and abiotic environment drives variability in levels of potential supply. In addition, other factors not considered here (both management-related and biotic or abiotic factors) likely also influence the potential supply of ecosystem services. At the forest stand scale, guidelines that aim to increase multifunctionality need to consider the existing trade-offs between ecosystem services, in order to set realistic goals. This is especially true for forest types aiming at timber production, as productivity showed synergies only with a small number of ecosystem services. Nevertheless, multifunctionality within forest types can possibly be improved via changes to those forest structural attributes that have strong effects on the potential supply of particular ecosystem services.

High levels of multifunctionality can also be achieved at the landscape scale by promoting those forest types that provide complementary subsets of ecosystem services. Within such an optimized landscape, each forest type would provide high levels of a subset of the desired ecosystem services at the stand level, theoretically resulting in high levels of all ecosystem services at the landscape scale. A mix of forest types will also increase the chances that one or more forest types provide additional ecosystem services, which could not be assessed or predicted in this study (e.g. scenic beauty, recreational value or hunting). Such a large-scale approach needs to take into account both the possibilities and limitations to establishing a particular forest type in a region, as well as the spatial variability of demands for ecosystem services (Burkhard et al. 2012). Hence, achieving multifunctional forest landscapes which supply ecosystem services where they are needed, and which consider both managers and users of ecosystem services, requires long-term and spatially explicit planning of forest management over large areas. The diversity of forest ownership within a given area complicates such top-down planning approaches but may in turn already provide multifunctional forest landscapes where differently managed forest stands provide different ecosystem services.

Conclusions

Both the assessment and the management of ecosystem services require a consistent mapping of forest structural attributes, across scales. Forest inventories are an excellent source of information for forest management on a nationwide scale and can be used to facilitate reporting of ecosystem functions and potential supply of ecosystem services with the help of predictive models such as ours. The lack of empirical support for deriving ecosystem service maps from land cover/land use maps has for some time been highlighted as a pressing challenge for policy making (Eigenbrod et al. 2010; Meyer et al. 2015). In contrast to maps based on expert knowledge or coarse land cover classes, our models provide a new empirical and semi-mechanistic basis for ecosystem service supply maps (Naidoo et al. 2008; Haines-Young et al. 2012). Hence, we encourage scientists and stakeholders to test our models for additional forest types and regions where the same or similar forest structural and compositional attributes are available and to complement them with existing predictive models for additional ecosystem services. Using forest inventories to assess and monitor the supply of ecosystem services alongside timber stocks and productivity can help guide a long-term and sustainable shift from a pure timber-focused to a truly multi-purpose forest management. The approach used here may also support much-needed payment for ecosystem service schemes to reward forest owners for providing multiple non-timber services.