Supplementary material and methods
DATA COLLECTION
In each forest patch, in both surveys (2003-2005 and 2014), we surveyed fungi and deadwood both in a survey plot and in the remaining patch area. The survey plot was of the size 20 m x 100 m and subdivided into survey cells of 20 m x 20 m (Fig. 1 in the main text). Within the cells we inventoried all deadwood with a diameter at breast height or base (depending on the deadwood type; referred to henceforth simply as diameter) ≥5 cm and length ≥1.3 m. For each spruce deadwood resource unit, we recorded the presence of fruit bodies of the focal species, spruce log diameter, length and decay stage. In the rest of the forest patch, i.e. in the remaining area outside the plot, we followed the same procedure but used a log diameter threshold of ≥15 cm. The lower resource-unit resolution in the remaining patch area was for covering the whole patch and for reaching a greater total number of patches than would have been possible if small-diameter deadwood was included everywhere.
All patches were dominated by Norway spruce (Picea abies ) and covered a range of forest types: clear-cuts with retention trees (53 patches, 16 of which had a plot, 16 x 5 cells, Fig. 1 in the main text), age of dominant trees 5-20 years at the time of the second survey; brook-side and herb-rich woodland key habitats (56 patches, 56 plots, 56 x 5 cells), age 64-163 years, representing natural-like set-asides; and managed forests (65 patches, 65 plots, 65 x 5 cells) of varying age, 66-134 years, and intensity of management. Patch area ranged from 0.2 ha to 4.3 ha (mean 1.1 ha). We recorded the forest stand variables living spruce volume and stand age in the first survey. Stand age in the second survey was the number of years in between the two surveys added to the age in the first survey. Volume of living stand in 2014 was predicted based on the stand characteristics measured during the first inventory using the MOTTI forest stand simulator (Salminen et al. 2005, Hynynen et al. 2014). MOTTI is a stand-level forest management and decision support tool which includes tree-level models for predicting growth and mortality. Prediction models for stand dynamics and underlying data sets are presented in detail in Hynynen et al. (2014).
Focal polypore species
Our two focal polypore species have spruce deadwood as their main resource (Berglund et al. 2011), are very rarely found in patches 20-60 years old (because such forests lack suitable deadwood; H. Moor et al., unpublished data) which our empirical data did not cover; are absent on stumps and living trees (Berglund et al. 2011), and were expected to have differing occupancies and colonisation-extinction rates because of differing ecologies (see below). Phellinus viticola has fairly large and easily recognized perennial fruit bodies leading to high detectability (Halme and Kotiaho 2012). Phellinus ferrugineofuscus has annual-biennial fruit bodies but as they are large and also well recognizable long (up to years) after the fruit bodies have died, this species also has high detectability. Both species further usually produce many fruit bodies per log during many years so that the total fruiting period is several years, even >10 years (Ovaskainen et al. 2013). Phellinus viticola occurs on both small- and large-diameter logs (Juutilainen et al. 2011, Nordén et al. 2013). Phellinus ferrugineofuscus occurs more frequently on larger logs but occasionally also on the thinner logs (Nordén et al. 2013). Phellinus ferrugineofuscus is ecologically relatively specialised while P. viticola is rather generalised in its resource use. This likely explains the different proportions of occupied patches by the two species – P. ferrugineofuscus in 17% andP. viticola in 24% of the patches in 2014.
COVARIATES EXPLAINING OCCUPANCIES AND METAPOPULATION DYNAMICS
We hypothesised that covariates collected for different spatial modelling scales would explain the occupancies and colonization-extinction dynamics of the focal species. The resource unit covariates were density (ha-1) and volume (m3 ha-1) of spruce logs, their mean diameter and stage of decomposition. The patch-level covariates were living spruce volume, age of the living stand (range 4-206 years), wetness, slope and aspect, measured as in Mair et al. (2017). Landscape-scale covariates were connectivity to potential dispersal sources in the surrounding landscape, and annual temperature and precipitation (Mair et al. 2017). For our relatively rare and specialized focal species we hypothesised positive effects of all resource-, patch- and large-scale covariates on occupancy and colonization probability but negative effects on extinction probability. Two exceptions are stage of decomposition and wetness for which we expected a unimodal relationship. These hypotheses are supported by a large number of studies (e.g. Heilmann-Clausen and Boddy 2008, Abrego et al. 2017, Nordén et al. 2018). Connectivity was calculated using satellite-based estimates of spruce volume and forest age in the form of 25 m grid cell layers (Tomppo et al. 2008). We followed the procedure of Mair et al. (2017), weighing the volumes of spruce in surrounding grid cells up to 20 km from the focal patch by\(\mathrm{\exp}\left(-\alpha d\right)\), where \(d\) is the distance to the focal patch and alpha is the relevant spatial scale of connectivity. We tested two different thresholds in forest age for inclusion of grid cells in the connectivity calculation (grid cells with stand ages ≥100 years or ≥120 years). We also tested replacing the spruce volumes in the connectivity measure with the presence or absence of forest exceeding the age threshold in each grid cell. As there is limited information available on species’ dispersal distances we further tested three different values for \(\alpha\) (0.1, 0.2 and 1), representing a mean dispersal distance of 10, 5 and 1 km, respectively. All covariates were log transformed, centred and standardised to have a mean of zero and a variance of unity.
Estimating the detection probability
One potential pitfall in models based on survey data is failing to account for unobserved species presences which can result in biased estimates of population trends and extinction probabilities (Kéry et al. 2006, Rhodes et al. 2006). Non-detection has likely a minimal effect on our projections as we accounted for it in several ways. Imperfect detection because of species characteristics, surveyor behaviour or sampling design can result in incorrect conclusions about population trends if not accounted for with the sampling and model design (MacKenzie et al. 2003, Altwegg et al. 2008). Firstly, in fungi, the species might not be fruiting at the time of the survey even if it is present as mycelia inside the dead tree. However, our focal species have similar mycelial and fruit-body prevalences as well as short time delays from mycelial colonisation to fruiting (Ovaskainen et al. 2013) which makes fruit body data meaningful and reliable. Secondly, the species might go non-detected even if it is present and fruiting, but we accounted for this in our models. Finally, a subset of the resource units may be included in the survey (e.g. based on minimum diameter), which we show can result in an opposite conclusion as compared with including all resource units: P. viticola was predicted to decrease in production forests when its smallest resource units were excluded, but predicted to increase when all the resource units were included.
A further hierarchical level in Bayesian modelling allows the unobserved “true” state of the population to be linked to the observed state via a model for the observation process and thus accounting for detectability (Harrison et al. 2011). For estimating the detection probability of the fruit bodies of the two polypore species, a larger survey plot of 4 ha (64 cells of the size 25 m x 25 m) in an old-growth forest patch in southern Finland was first surveyed in 2002 and partly (15 cells) resurveyed twice in 2014. The principle for detectability data collection is similar as in the general survey, and hence the detectability estimate obtained was used in the modelling, see below. In 2014, a field team of two experts first carried out a regular resurvey without knowledge about the subsequent second resurvey. In the second resurvey, the ‘intensive resurvey’, the task was to do as thorough a survey as possible, taking the time needed, in a team of four experts that now knew the purpose of the intensive resurvey was to collect data for estimating detectability. We fit detectability models to these data assuming the occupancies recorded during the intensive resurvey represented the true occupancies of the polypore species encountered during the surveys. As the density of deadwood is considerably larger overall in this old-growth forest than in a typical managed forest, we included deadwood density as a covariate in the detectability models. The mean detection probability across all polypore species recorded in the intensive resurvey was estimated to be 0.9 for the mean deadwood density in managed forests, and >0.9 in forests with a low density of dead wood. In the colonisation-extinction and occupancy models, we therefore used a patch-specific detection probability 0.9.
MODELLING OCCUPANCY AND COLONISATION-EXTINCTION
For each species, we fitted hierarchical Bayesian state-space models to the presence-absence data of the species at three spatial modelling scales (cell, plot and patch) and two deadwood resource resolutions (diameter ≥5 cm or ≥15 cm). We included covariates collected for different spatial modelling scales that we hypothesised would explain the occupancies and colonization-extinction dynamics of the focal species. For the comparison with the colonisation-extinction models we also fitted occupancy models. We used data from the second survey (174 patches surveyed both in 2003-2005 and 2014), and, for obtaining more reliable estimates of occupancy (but not colonization-extinction requiring two surveys) on clear-cuts, especially in relation to clear-cut age, we further included data from 19 patches from the first survey (2003-2005) that were not re-surveyed in 2014, and 22 patches that were cut between the two surveys and therefore only surveyed in the first survey.
In the modelling of species occupancy in clear-cuts, we only tested the effects of stand age and connectivity, and used all first-survey and resurvey data from clear-cut patches. We hypothesized (i) a negative age effect with time since cutting, and (ii) a positive effect of connectivity to surrounding dispersal sources. We did not test for the effect of deadwood quantity in clear-cut patches as most of the dead wood in these patches is made up of logging residues which, based on our data, are very seldom used by the focal species.
The response variables in the models (probabilities of occurrence, colonisation and extinction) were adjusted to a time scale of 10 years (only colonisation and extinction concerning change through time) and a spatial modelling scale of one hectare (Appendix S1). This was conducted by offsetting for the differing numbers of years between the surveys and the varying survey areas (cell, plot or patch). In the cell-level models, we included a random effect for the patches, which was assumed to be normally distributed with mean zero and a variance that was estimated during model fitting. In the projections we assumed a colonisation probability of zero and extinction probability of unity at such locations.
We fitted separate models at each of the cell-, plot- and patch scales and we detail the cell-scale model. The plot- and patch-scale models are simplifications of this model with the finest level indexing instead at the plot and patch scale, respectively.
For the cell-scale models we define \(Z_{i,j,t}\) as the occupancy state of the focal species in cell \(i\) in patch \(j\) during survey period\(t\). We assume that\(Z_{i,j,t}\ \sim\ \text{Bernoulli}\left(\psi_{i,j,t}\right)\). For the first survey period \(\psi_{i,j,t}=\omega\) and in the following period:
\(\psi_{i,j,t}=\left(1-Z_{i,j,t-1}\right)c_{i,j,t}^{*}+Z_{i,j,t-1}\left(1-e_{i,j,t}^{*}\right)\),
where \(c_{i,j,t}^{*}\) and \(e_{i,j,t}^{*}\) give the colonisation and extinction probabilities, respectively. They have been offset to correct for the different numbers of years between the surveys and the different cell areas (20 m x 20 m or 25 m x 25 m), such that:
\begin{equation} c_{i,j,t}^{*}=1-\left(1-c_{i,j,t}\right)^{n_{i,j,t}a_{i,j,t}}\nonumber \\ \end{equation}\begin{equation} e_{i,j,t}^{*}=1-\left(1-e_{i,j,t}\right)^{\frac{n_{i,j,t}}{a_{i,j,t}}}\nonumber \\ \end{equation}
where \(n_{i,j,t}\) gives the number of years between the surveys divided by 10 (i.e. scaled by the typical number of years) and\(a_{i,j,t}\) gives the cell area divided by 400 (i.e. scaled by the typical cell size in m2). For a cell with ten years between the two surveys and an area of 400 m2,\(c_{i,j,t}^{*}=c_{i,j,t}\) and \(e_{i,j,t}^{*}=e_{i,j,t}\). If the forest in the cell had not been clear-cut we assumed that:
\begin{equation} \text{cloglog}\left(c_{i,j,t}\right)=\delta_{2}+\sum_{k}{\beta_{k}^{(C)}X_{k,i,j,t}^{(C)}}+\sum_{l}{\beta_{l}^{(P)}X_{l,j,t}^{(P)}}+\gamma_{j,t}\nonumber \\ \end{equation}
and that \(\text{cloglog}\left(e_{i,j,t}\right)=\varepsilon_{2}\). We chose to use the complementary log-log link function,cloglog, as due to its asymmetrical nature it is better suited than the more conventional logistic link function to cases where the probabilities are very large or very small (Golding et al., in preparation). Due to data limitations we could not include covariates or random effects in the models for the extinction probabilities or the colonisation probability for cells in forest patches that had been clear-cut (either before the first survey or between the two survey events) and intercept only models were used in these cases. Specifically, we estimated separate parameters for clear cuts,\(\mathrm{\text{cloglog}}\left(c_{i,j,t}\right)=\delta_{1}\) and\(\text{cloglog}\left(e_{i,j,t}\right)=\varepsilon_{1}\). The \(k\)cell-scale covariates used in the model for the colonisation probability are given by \(X_{k,i,j,t}^{(C)}\) with corresponding parameter values\(\beta_{k}^{(C)}\). The \(l\) patch-scale covariates are given by\(X_{l,i,j,t}^{(P)}\) with corresponding parameter values\(\beta_{k}^{(P)}\). We also include a patch-scale random\(\gamma_{j,t}\ \)which was assumed to be normally distributed with mean zero and variance \(\sigma^{2}\). In cases where colonisation and extinction probabilities were observed to be 0 and 1, respectively (see Table 1), we used these probabilities in the simulations.
Finally, we define \(Y_{i,j,t}\) as the observed occupancy state of cell\(i\) in patch \(j\) during survey period \(t\). For the observation model we assume that\(Y_{i,j,t}\ \sim\ \text{Bernoulli}\left(Z_{i,j,t}p\right)\) where\(p\) gives the detection probability. This detection probability was estimated based on an intensive control study (see Estimating the detection probability above) and a value of 0.9 was used.
For the cell-scale models it was necessary, due to statistical non-independence, to include a random patch effect term. At the higher spatial modelling scales (plot and patch) it was not possible to include such random effects. This random effect term has the advantage that it will capture between patch variation that is not explained by the covariates included in the models, for instance due to unmeasured covariates, habitat differences at the time of establishment but not necessarily at the time of the inventories, or competitive exclusion by species already established on some patches (Berglund et al. 2009b). True absences on some suitable patches can also be a consequence of stochasticity in the dynamics of either the species (Snäll et al. 2005) or their host trees (Jonsson 2000) or both. Not accounting for such variability has been shown to bias model-based predictions (Berglund et al. 2009a), especially for resource-unit restricted species such as wood-decomposer PROJECTING fungi.PROJECTING FOREST CONDITIONS THE COMING CENTURYTHE COMING CENTURY
We made projections of the distribution of the resource of our focal species on 17,383 National Forest Inventory (NFI) plots spread across the boreal zone of Sweden between 2010 and 2110 under a ‘business as usual’ scenario utilizing the recent nationwide Forestry Scenario Analysis by the Swedish Forest Agency (Claesson et al. 2015, Eriksson et al. 2015). In this scenario, 16% of forest land is protected in state reserves, patches voluntarily protected by individual landowners or areas set aside from forestry, retention patches. We refer to the collection of these non-production patches as set-asides in what follows. The set-asides are situated in forest landscapes otherwise dominated by clear-cutting forestry. The NFI plots are circular and have a diameter of 7 or 10 m but for all the projections we first scaled all variables to values at one hectare.
Projections of forest management and dynamics were made using the multiple-use forestry planning tool Heureka (Wikström et al. 2011). For the projections the RegWise software (version 2.2) was used. RegWise is a software component within the Heureka suite of forest decision support tools (Wikström et al. 2011). The core of Heureka is made up of empirical individual models for tree growth (see Fahlvik et al. 2014), ingrowth (Wikberg 2004) and mortality (Elfving 2014) together simulating tree layer development in five-year time steps. In RegWise forest management actions are steered by a rule-based simulation framework. The harvest level – and consequently the development of the forest – is controlled by specified management programmes (silviculture and harvest activities) and the actual growth in each time period. Different management programmes are specified for, e.g., different forest owner categories, tree species mixtures and site factors. Tree retention practices at final felling were also included. The forest projection output data included stand age, a covariate used in the modelling, living spruce volume which was utilised in our connectivity calculations and spruce deadwood production generated by natural tree mortality. We initialized the projections with the amounts and properties of spruce deadwood observed on the NFI plots during 2008-2012. We simulated decomposition of the deadwood based on the one-time chronosequence method described in Harmon et al. (2000). When deadwood enters the fourth stage of decomposition (based on the Swedish decay stage classification (Sandström et al. 2007), it is removed from the simulations as our focal species do not produce fruit bodies on such highly decomposed wood. We do not differentiate between standing (snags) and downed (logs) deadwood in the simulations but used the proportion of spruce logs (out of all dead spruce) in the NFI data to determine how much to multiply the deadwood quantities by in the covariates for our focal species. These multipliers were estimated as 0.59 for the volume of spruce dead wood and 0.64 for the density of dead spruce trees. Forest connectivity projections were done as in Mair et al. (2017).
All the model parameters were given priors with a mean of zero and a standard deviation of 0.01, except for sigma which was given a uniform prior between zero and five. In all cases three MCMC chains were run for 3000 iterations with a burn-in of 1000 iterations and a thinning of 100. Convergence was assessed visually based on the results of the three chains. In all cases convergence was reached after 3000 iterations.
PROJECTING POLYPORE OCCUPANCIES THE COMING CENTURY
We projected the occupancy dynamics of our focal species for 2010-2110. For each polypore species the final fitted occupancy model was utilised to initialise the occupancy states in the first time step, here 2010. As the average time between our surveys was 10 years, we used 10-year time steps to simulate the subsequent colonisation and extinction dynamics in the NFI plots until 2110 using the final fitted colonisation-extinction model with its estimated parameter values. Offsets were used on the occupancy and colonisation and extinction probabilities to account for the different sizes of the NFI plots (7 and 10 m radii). Finally, area factor scaling, accounting for different densities of NFI tracts across Sweden, was used to compute means and totals of these hectare-scale values. Although our model fitting, projections and decomposition modelling are in terms of deadwood with diameter ≥5 cm at the finest resolution, the NFI data used for initialising the projections are in terms of deadwood with diameter ≥10 cm. Visual inspection of the modelling output showed that after ten years the effect of the missing small pieces (<10 cm) on the deadwood projections was minimal and we thus show the results between 2020 and 2110.