Stand age and climate influence forest ecosystem service delivery and multifunctionality

We examine how levels of multiple ecosystem services (ESs) change with succession in forests with different tree species composition. More specifically we ask how ecosystem age interacts with environmental conditions to regulate ES delivery. Using the nationwide Swedish forest inventory, comprising boreal and temperate regions, we investigated how levels of six provisioning, regulating, recreational, and/or cultural forest ESs changed with forest age (10–185 years) in stands of different tree species composition. We also tested whether the number of ESs delivered (i.e. multifunctionality) changed substantially with stand age, using different threshold levels for ES delivery. Accounting for environmental conditions and stand properties, we found that levels of single ESs changed with stand age. Tree biomass production usually peaked in young to medium aged stands. In contrast, production of berries and game, and services related to biodiversity, were typically highest in old stands (120–185 years). Consistent with this strong temporal tradeoff, multifunctionality at lower threshold levels increased with stand age in most monocultures and mixtures, with the highest multifunctionality being reached somewhere between 100 and 185 years, depending on tree species composition. This was not evident for the highest threshold ES level (the top-20%), however. Moreover, multifunctionality usually decreased with warmer climatic conditions, with the exception of spruce–pine–birch mixtures. Taken together, our results show that a reduced forest age, e.g. due to forestry targeting early harvest of stands, most likely would limit the delivery of several ESs valued by society and result in less multifunctional forests. To maintain the capacity of forests to deliver high levels of multiple ESs, the role of stand age and tree species composition should be considered in decisions on how to manage future forests.


Introduction
The number and levels of ecosystem services (ESs) delivered by a forest depend on several stand properties, such as tree species richness (Gamfeldt et al 2013), species composition (Gamfeldt et al 2013, Baeten et al 2019, Jonsson et al 2019, relative abundance (Jonsson et al 2019), structural complexity (Felipe-Lucia et al 2018, Strengbom et al 2018, and environmental conditions, such as light levels and soil chemistry (Gamfeldt et al 2013, Ratcliffe et al 2017, Felipe-Lucia et al 2018, Strengbom et al 2018. All of these forest properties, and forest floor plant community composition, change rather predictably with forest succession (Wardle et al 2012, von Arx et al 2013, Kumordzi et al 2015. Earlysuccessional forests are usually single-layered, especially after clearcutting, with smaller trees at high stem density. As such stands grow taller, the initially light ground-floor conditions get increasingly darker. Mid-successional forests are more multilayered, with closed canopies but more open ground level, allowing for shade-tolerant species to exist. Late-successional forests tend to be multi-layered, but also exhibit gap dynamics that reduce stem density and increase light and habitat heterogeneity at the ground level (Mccarthy 2001). As canopy trees are strong competitors compared to many groundfloor plant species, the described successional changes should strongly affect ES delivery in forests, especially services depending on a rich and dense ground-floor layer (e.g. berry production, understory plant species richness, mushroom production, food for game) (Hedwall et al 2019a ).
Despite the fact that successional changes in forest properties are rather well understood, the effects of succession, as indicated by stand age, on delivery of multiple ESs have rarely been quantified (but see Sutherland et al 2016, Jönsson andSnäll 2020). In managed forest landscapes, stand age is a key variable that, in addition to tree species composition, is determined by the manager. In boreal forests, natural stands can be several hundred years old, whereas managed forests rarely are allowed to get older than 100 years, and often only 70-80 years, before being harvested (Esseen et al 1997, Larsson andDanell 2001). If maximum delivery of some ESs occurs at later stand ages, management that restricts stand age in the landscape will have large-scale consequences for the delivery of different ESs, and thus likely for forest multifunctionality. Hence, understanding how the delivery of different forest ESs change as stands age, and if these changes differ among stands containing different tree species, is important in the formulation of policies for sustainable forest management (Díaz et al 2015, Sutherland et al 2016, Jonsson et al 2019, Jönsson and Snäll 2020. Here we use data from a nationwide forest inventory and a space-for-time substitution approach, to examine how the delivery of six ESs, that are classified as provisioning, regulating, recreational, and/or cultural (Jonsson et al 2019), changes with stand age (10-185 years) in stands of different tree species composition. We also investigate if the number of ESs delivered (i.e. multifunctionality) changes with stand age. Our main hypotheses are that (1) both the number and levels of ESs change with stand age, (2) such changes differ depending on tree species compositions, and (3) there are temporal tradeoffs between tree biomass production and other forest ESs, implying that all services cannot be maximized at the same time.

Materials and methods
We used data from the Swedish National Forest Inventory (NFI) and the Swedish Survey of Forest Soils and Vegetation (SFSV; www.slu.se/markinventeringen).
For detailed information about these data and our criteria for data selection, see supplementary information.

Ecosystem services
Boreal and temperate forests are extensive biomes that globally provide a wide range of ESs to humans (Hansen et al 2010, Gauthier et al 2015. Tree biomass production was estimated as the yearly change in tree biomass (kg m −2 yr −1 ), calculated over a period of 5 years, and for all tree individuals with a diameter ⩾40 mm at breast height (1.3 m). For plots visited in the years 1999-2002, the baseline for biomass production was thus the years 1994-1997. Sample size for tree production was 4335 plots, and due to fallen trees, observed values could be negative (as was the case in 5 plots). For more detailed information, see supplementary information.
Soil carbon (C) storage was measured as the amount of carbon (g m −2 ) in the topsoil, which consisted of either purely organic horizons, e.g. mor layers (63%) or peat layers (21%), or less frequently of minerogenic A-horizons (16%). This is the part of the soil most affected by the current aboveground biota. To compensate for the conceptual difference in topsoil types, mean soil C stocks were set equal for purely organic soils (measured in the whole organic horizon down to a maximum depth of 30 cm) and minerogenic topsoils (measured in the top 10-cm horizon) (Gamfeldt et al 2013). The soil fraction <2 mm was analyzed. Soil sampling was carried out in around 50% of the inventory plots, totaling 1953 plots. Data were log transformed.
Bilberry production was measured as the percentage of each plot covered by bilberry, Vaccinium myrtillus. The cover of V. myrtillus is strongly correlated to annual production, and bilberry is one of the economically most important wild berry species in Northern Europe (Miina et al 2009). Bilberry sampling was carried out on around 50% of the inventory plots, totaling 2127 plots. Data were converted to proportions and arcsine transformed, after correcting for the area where berries could not grow, e.g. the area of stems and boulders.
Game production potential was calculated as the percentage cover of each plot that was covered by any of 15 plant species important for large herbivorous mammals, for example, moose (Alces alces) and roe deer (Capreolus capreolus) (Cederlund et al 1980). The cover of these field layer species was measured in the same way as the bilberry cover. Sample size was 2127 plots and data were log transformed after adding a constant of 0.01. For more detailed information, see supplementary information.
Understory plant species richness: The Swedish NFI includes about 270 ground vegetation species, and of these, we selected those with forest as main habitat (in total 141 species). We modeled the proportion of the richness of these species (i.e. number of species) in each plot in relation to the regional species richness of the same species, estimated as all of the 141 species occurring within a 600 × 600 km window centred on each plot (Gamfeldt et al 2013). Hence, this ES is a proportion between 0 and 1. Sample size was 4327 plots. For more detailed information, see Supplementary information.
We modeled the probability of dead wood occurrence using data on presence/absence of dead wood in each plot (minimum diameter for dead wood pieces to be counted equals 40 mm). We chose this variable, as the abundance of dead wood showed a highly skewed distribution in which small fragments were common and a few plots had very high abundances. Sample size was 4335 plots. For more detailed information, see supplementary information.
Number of ecosystem services (i.e. multifunctionality) was the number of services found present in each plot. To assess whether multifunctionality depends on the level of delivery, we used four different thresholds; 100% (where all values >0 for each service were included) and top-80%, top-50%, and top-20%, where these were classified as present (= 1) for each service, and the bottom-20%, bottom-50%, and bottom-80%, respectively, were classified as absent (i.e. zero).

Predictor variables
Temperature sum (or degree days; henceforth abbreviated into temperature) is a measure of local climate, and is in the NFI defined as the total accumulated mean daily temperature (in • C) over 5 • C during the vegetation period. Humidity was measured as the difference in mm between the addition of water via precipitation and the loss of water due to transpiration and evaporation during the vegetation period. For nitrogen (N) deposition, spatially distributed N deposition data in kg N ha −1 yr −1 (wet and dry) were downloaded from the Swedish Meteorological and Hydrological Institute website (www.smhi.se). Temperature, humidity, and nitrogen deposition are 'modeled' variables with a tract-scale resolution. pH of the topsoil was measured in water suspension. The C:N ratio was the ratio of total carbon to total nitrogen in the topsoil obtained from the Soil Inventory. Soil moisture was measured on a categorical scale of 1 (driest) to 5 (wettest). These categories correspond to a continuous scale of % volumetric soil water contents ranging from 0.15 to 0.35 (Warfvinge and Sverdrup 1995). Peat was modeled as peat (1) or non-peat (0) soil. Soil is classified as peat when no mineral soil is found within a 30-cm depth from the soil surface.
In this study, we focused on monocultures of Norway spruce (Picea abies), Scots pine (Pinus sylvestris), and birch (Betula spp.), and mixtures of these three species, and additionally 'other broadleaved' species. Spruce, pine, and birch are considered commercially viable in Sweden, and make up approximately 41%, 39%, and 12%, respectively, of the productive forest area in Sweden (Forest Statistics 2018). 'Other broadleaved' in our NFI subset dataset contains at least two, but not more than five species, and is mostly mixtures of oak (Quercus robur), aspen (Populus tremula), beech (Fagus sylvatica), and alder (Alnus spp.). Consequently, we ended up with the six stand types: spruce, pine, spruce-pine, spruce-birch, spruce-pine-birch, and spruce-pine-birch-other broadleaved. Birch monocultures and pine-birch mixtures were excluded from subsequent statistical analyses, due to poor representation across the entire stand-age gradient. Forest stand age was estimated in each plot as the basal area weighted mean age based on tree-ring counts from cores of at least three trees, independent on species identity. In the statistical analyses, and when plotting data/results, only stands with an age ⩽185 years were included, because too few data points beyond that age limit were present in our data. Stem density was the sum of two stem densities (i.e. counts of large trees and seedlings) in each plot. The number of stems taller than breast height (1.3 m) and with a diameter >40 mm were counted in each plot (radius of 10 m) and then divided by total plot area (314 m 2 ). This stem density was summed with the number of stems taller than breast height (1.3 m) and with a diameter <40 mm that were counted in a subplot (3.5 m radius) and then divided by the subplot area (38.5 m 2 ). As the distribution of stem density was right-skewed with a long tail, we used log transformed values in the statistical analyses.

Statistical analyses
To investigate temporal variation in ES delivery, we used a space-for-time substitution approach (Pickett 1989), and modeled each of the six ESs (y) as a function of stand age (x) and the other environmental variables, for the six different stand types. The full initial model used for every ES and tree species composition was: Ecosystem service ∼ stand age + stand age 2 + pH +pH 2 + N deposition + C : N ratio + soil moisture +soil moisture 2 + stem density + humidity +temperature sum + peat soil + stand age ×temperature sum + tract as random effect As total standing tree biomass is a forest characteristic that changes predictably with stand age, and as effects of age is our primary question, we did not include tree biomass as a covariate in the model.
The partial model y ∼ x + x 2 , including 95% confidence bands, of each final model, was plotted across the gradient in stand age as a proportion of the maximum observed delivery for each service, independently of species composition. The results for four of these stand types were also plotted in petal diagrams for the stand ages 30, 70, 110, and 150 years, using 'coord_polar' in ggplot (R Core Team 2019 ), for better visualization. These ages were chosen as they roughly represent the stand ages for pre-commercial thinning, lower limit for clearcutting, upper limit for clearcutting, and old-growth forest, respectively.
For each of the six stand types, we also modeled the number of ESs (y) as a function of stand age (x) and other environmental variables. However, as this subset dataset included far fewer plots than the set for each individual service, we limited the number of covariates to those most commonly included in the bestfitting models for the individual services. Hence, the full initial model for number of ESs was: Number of ecosystem services ∼ stand age + stand age 2 +pH + pH 2 + soil moisture + soil moisture 2 +stem density + temperature sum + stand age ×temperature sum + tract as random effect Model fitting was done for all ES threshold levels (see Number of ecosystem services above), to assess the sensitivity of using different levels of delivery to define multifunctionality. The partial model y ∼ x + x 2 , including 95% confidence bands, of each final model that included stand age as significant was plotted. However, as temperature turned out to be a predictor variable of major importance, we also plotted the partial model of each final model that included temperature as significant. There were also two cases where the interaction stand age × temperature was significant.
For the analyses, missing values (NAs) were removed, and for the models on individual services we used Akaike information criterion (AIC) to obtain the best-fitting model. Generalized linear mixedeffects model (GLME) was used to fit individual services ∼ stand age, using 'lmer' or (for probability of dead wood occurrence) 'glmer' in the R package lme4 (Bates et al 2015). GLME was used also to fit number of ESs ∼ stand age, but here we used 'glmmPQL' for data with negative binomial distribution in the R package MASS (Venables and Ripley 2002). As glmmPQL uses quasi-likelihood (rather than log-likelihood) functions, AIC is not available for model selection. Hence, the model selection was based on estimated standard errors and P-values for the coefficients associated with the variables. To plot the partial model for each individual service, we used stat_smooth, function 'lm' , in the R package ggplot2 (2.2.1.9000), with 95% confidence bands (R Core Team 2019). Number of plots per stand type and individual ESs are presented in tables S1-S6 (available online at https://stacks.iop.org/ERL/15/0940a8/ mmedia), and the same information for number of ESs is presented in tables S7-S12.

Results
Stand age had a significant effect on the level of ES delivery in 28 of the 36 investigated cases (figure 1, tables S1-S6). Tree biomass production often showed a hump-shaped relationship with increasing stand age, whereas bilberry production and probability of dead wood occurrence always increased with stand age, and were highest at stand ages above 150 years (figure 1). Game production potential was relatively high in young stands and often the lowest at 70-90 years. However, later it increased in spruce monocultures and in spruce-birch mixtures, while it decreased in pine monocultures ( figure 1).
For all tree species compositions, the highest levels of tree biomass production were found in stands <80 years, when the other ESs-especially bilberry production and dead wood occurrence, but also game production potential-mostly showed low levels (figures 1 and S1). This early increase in tree biomass production (often followed by a peak in soil C storage) was most apparent for spruce monocultures and spruce-birch and spruce-pine-birch-other mixtures, up to a stand age of 50-80 years of age, after which a decrease in tree biomass production followed. Among these stand types, the highest tree biomass production was found in spruce-birch mixtures, closely followed by the other two stand types and the youngest spruce-pine stands (figure 1). In stands above 150 years of tree age, tree biomass production was close to zero, independently of tree species composition (figure 1).
Besides stand age, a range of environmental conditions was important for levels of the six ESs investigated. Although there were differences among stand types and services in terms of which variables were included in the final models, some general patterns emerged (tables S1-S6), with temperature most frequently included (28 models), followed by soil moisture (20 models), and pH (19 models). Not surprisingly, for tree biomass production, temperature, soil moisture, and stem density (all positive) were most frequently included in the models, whereas temperature and pH (both negative) were most frequent for bilberry production.
For soil C storage, all models showed a positive relationship to soil moisture and peat soil, and temperature (positive) was included in 5 out 6 models (tables S1-S6). In contrast, for game production potential, there were negative relationships to temperature (4 out of 6 models) and N deposition, but positive relationships to stem density (3 out of 6 models). Independently on stand type, increasing pH always promoted understory vegetation species richness, though increasing soil moisture and decreasing temperature were also important. Lastly, the probability of dead wood occurrence decreased with temperature and N deposition in 5 out of 6 best-fitting models (tables S1-S6).
Multifunctionality increased with stand age in 8 out of 24 possible cases, and marginally significantly so in one additional case (tables S7-S12), reflected in the gradually increasing number of ESs delivered with increasing stand age (figure 2). Not surprisingly, multifunctionality also decreased with Figure 1. Level of ecosystem service delivered, i.e. proportion of maximum delivery dependent on tree species composition and stand age, for tree biomass production, bilberry production, top-soil C storage, game production potential, understory vegetation species richness, and dead wood occurrence, across a stand age gradient for different selected tree species compositions. Plotted lines are partial regression models with 95% confidence bands (shaded area) from final GLME models, accounting for other explanatory variables (supplementary tables S1-S6). Stand age had no significant effect on top-soil C storage in spruce, spruce-birch, and spruce-pine-birch-other stands, on understory vegetation species richness in spruce, pine, and spruce-pine-birch-other stands, or on game production potential in spruce-pine and spruce-pine-birch-other stands.
increasing threshold, and no significant relationships between stand age and number of ESs were found at the highest threshold (i.e. top-20%), indicative of increasing tradeoffs among services with increasing level of delivery. Multifunctionality in spruce monocultures showed the strongest relationship to stand age. The mean number of ESs delivered increased by about two along the stand-age gradient, for both the top-80% and the top-50% thresholds (figure 2). With the exception of tree biomass production and soil C storage, individual ESs (tables S7-S12) and multifunctionality decreased with temperature (i.e. in warming climates) along the latitudinal gradient studied ( figure 3). However, in spruce-pine-birch, multifunctionality increased with increasing stand age at high temperature, i.e. in warmer climates (figure S2), as revealed by a significant stand agetemperature interaction (tables S11a and S11c).

Discussion
We found that there was a general increase in multifunctionality with increasing stand age in the most common tree species compositions in Scandinavian boreal and temperate forests. This multifunctionality increased northwards, to forests in colder climates. Most ESs peaked in forest stands older than 120 years, with the exception of tree biomass production which was highest at stand ages of 50-70 years. It has been shown previously that forest attributes, such as overall tree species richness or specific tree species composition (e.g. Gamfeldt et al 2013, Felipe-Lucia et al 2018, Baeten et al 2019, Jonsson et al 2019, specific forest management methods (Strengbom et al 2018), or whether a forest is managed or not (Peura et al 2018), determine the delivery of multiple forest ESs. Our study adds to these insights by showing that stand age is at least as important as certain stand characteristics in determining not only the levels of individual ESs, but also the number of services delivered (i.e. multifunctionality) (see also Sutherland et al 2016, Jönsson andSnäll 2020). This is highly relevant information for sustainable forest management, because restricting stand ages by clearcutting medium-aged and old stands is one of the strongest impacts of modern forestry.
Our results show that in most stand types, one or several ESs are maximized only after the typical age of clearcutting in boreal forests. Specifically, in boreonemoral Scandinavia, the rotation length, i.e. the time between clearcutting events, is usually between 70 and 100 years, depending on climate, site productivity, and the focal tree species (i.e. usually Scots pine or Norway spruce). Since most of the studied ESs were at their maximum after 120 years, the decreased stand age caused by modern largescale forestry has inadvertently resulted in forest landscapes that have lost much of their multifunctionality, as suggested in a simulation study by Triviño et al (2017). Reduced rotation lengths, which have sometimes been suggested as measures to increase harvest volumes, would therefore be expected to further reduce forest multifunctionality. Conversely, an extended rotation length would increase the levels of several ESs, though at the expense of reduced tree biomass production. As all forest services per definition are important for human well-being, it is important to investigate ways as to how management can mitigate such tradeoffs (Howe et al 2014). One way could be to manage forests to contain certain tree species mixtures at specific relative abundances (Jonsson et al 2019), but our results imply that levels of ESs provided by these mixtures still would change as the forest ages. This means that to maintain high delivery of multiple ESs s on the landscape scale, and forests that are locally highly multifunctional, sufficiently large and numerous areas of old forests should be set aside, in addition to increasing the rotation lengths in part of the managed forest landscape.
Resource competition between trees and other plant species is likely the reasons why old forests deliver higher levels of some ESs, and particularly those linked to cover and diversity of the understory vegetation (i.e. bilberry production, understory vegetation species richness, and game production potential; cf Widenfalk and Weslien 2009). Especially younger (<80 years) undisturbed spruce forests are highly dense and therefore outcompete understory vegetation for light (e.g. Strengbom et al 2018, Hedwall et al 2019a, Hedwall et al 2019a, 2019b. Accordingly, with increasing stand age beyond 80 years, the largest increase of some individual ESs (e.g. berry production), and of multifunctionality, took place in spruce monocultures and spruce-birch mixtures, likely because of gap dynamics in the old sprucedominated forests (Mccarthy 2001). However, in stands containing three or more tree species, there was no increase in understory plant species richness with stand age, indicating light limitation for understory vegetation also in older stands (Hedwall et al 2019a). On the other hand, in these stands, understory richness was relatively high in young forest despite comparatively high tree biomass production, suggesting that a tradeoff between wood production and understory plant richness in managed forests <80 years old can be mitigated in stands of more than two tree species. It is therefore possible that certain management actions, such as thinnings that increase the penetration of light to the ground, create conditions that favor some ESs (Bauhus et al 2009, Widenfalk and Weslien 2009, Strengbom et al 2018. Nevertheless, a production forest, albeit thinned, lacks several other old-growth characteristics, such as an uneven age structure, that are important for ES delivery (e.g. Felipe-Lucia et al 2018) but difficult to manage for (Bauhus et al 2009).
In Sweden, large grazing animals, and especially moose (Alces alces), are often considered a problem for forest management focusing on Scots pine, as the shoots on young pine are a preferred food resource. This is especially so in the absence of more palatable deciduous tree species (Bergqvist et al 2012, Milligan andKoricheva 2013). Accordingly, our results do indeed show that young pine monocultures have large game production potential, which gradually decreases as this stand type grows older. More surprisingly, however, also old spruce monocultures and spruce-birch mixtures showed high game production potential-almost 50% higher than young pine monocultures. This is likely explained by gap dynamics in the oldest spruce-dominated forests resulting in more diverse and abundant ground-and shrub-layer vegetation, which include preferred grazer food species. Figure 2. Number of ecosystem services delivered by forest stands, across a gradient of stand age, independent on level of delivery (100%), including only stands that delivered services in the top-80%, including only stands that delivered services in the top-50%, and including only stands that delivered in the top-20%. Only tree species compositions with statistically significant relationships between number of ecosystem services delivered and stand age are presented (dashed line: P < 0.1). Shaded area represents 95% confidence intervals. Plotted lines are partial regression models with 95% confidence bands (shaded area) from final GLME models, accounting for other explanatory variables (supplementary tables S7-S12).
Besides intensified forest management, climate change poses large, and largely unknown, challenges to forests and the services they deliver (Seidl et al 2017, Astrup et al 2018. This may be especially true for boreal forests, as the increase in temperature is projected to be the greatest at norther latitudes (IPCC 2014). In the light of this, our finding of decreasing multifunctionality with increasing temperature is highly relevant, as it suggests that future warming of boreal and northern temperate forests may lead to decreased multifunctionality. Such a decrease would likely be caused by increased density and height of trees in a warmer climate, and subsequent decreases in light availability on the forest floor, resulting in loss of light-demanding plant species (Strengbom et al 2018, Hedwall et al 2019bis, whose diversity associated with the delivery of several ESs. However, our results also suggest that management for attributes of mixed, old stands could mitigate such adverse effects of warming on forest multifunctionality. Thus, in addition to our earlier call for adopting management strategies that mitigate tradeoffs among services (Howe et al 2014), it is also important to promote forest characteristics that are better able to withstand adverse impacts of future climate regimes (Astrup et al 2018).
In summary, our results show that to meet current and future human demands on ecosystems to deliver a range of ecosystem services (Daily 1997, IPBES 2019, larger areas of old forest should be set aside, and forest attributes that promote the delivery of multiple services simultaneously (e.g. Felipe-Lucia et al 2018) should be managed for in production forests. A good candidate for such an attribute is Figure 3. Number of ecosystem services delivered by forest stands against temperature sum, a measure of local climate, independent on level of delivery (100%), including only stands that delivered services in the top-80%, including only stands that delivered services in the top-50%, and including only stands that delivered in the top-20%. Only tree species compositions with statistically significant relationships between number of ecosystem services delivered and temperature are presented (dashed line: P < 0.1). Shaded area represents 95% confidence intervals. Plotted lines are partial regression models with 95% confidence bands (shaded area) from final GLME models, accounting for other explanatory variables (supplementary tables S7-S12). canopy heterogeneity, either caused by natural gap dynamics in old stands or by the mixing of certain tree species in younger stands. However, exactly what the most important attributes are, and how to best mitigate the observed strong tradeoff between tree biomass production and most other services, remains questions to be resolved in future studies.