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.