DOI: 10.5061/dryad.2ngf1vhq5
title: Multiple drivers of large‐scale lichen decline in boreal forest canopies

description: Thin, hair-like lichens (Alectoria, Bryoria, Usnea) form conspicuous
 epiphyte communities across the boreal biome. These poikilohydric
 organisms provide important ecosystem functions and are useful indicators
 of global change. We analyse how environmental drivers influence changes
 in occurrence and length of these lichens on Norway spruce (Picea abies)
 over 10 years in managed forests in Sweden using data from &gt;6000
 trees. Alectoria and Usnea showed strong declines in southern-central
 regions, whereas Bryoria declined in northern regions. Overall, relative
 loss rates across the country ranged from 1.7% per year in Alectoria to
 0.5% in Bryoria. These losses contrasted with increased length of Bryoria
 and Usnea in some regions. Occurrence trajectories (extinction,
 colonization, presence, absence) on remeasured trees correlated best with
 temperature, rain, nitrogen deposition, and stand age in multinomial
 logistic regression models. Our analysis strongly suggests that industrial
 forestry, in combination with nitrogen, is the main driver of lichen
 declines. Logging of forests with long continuity of tree cover, short
 rotation cycles, substrate limitation and low light in dense forests are
 harmful for lichens. Nitrogen deposition has decreased but is apparently
 still sufficiently high to prevent recovery. Warming correlated with
 occurrence trajectories of Alectoria and Bryoria, likely by altering
 hydration regimes and increasing respiration during autumn/winter. The
 large-scale lichen decline on an important host has cascading effects on
 biodiversity and function of boreal forest canopies. Forest management
 must apply a broad spectrum of methods, including uneven-aged continuous
 cover forestry and retention of large patches, to secure the ecosystem
 functions of these important canopy components under future climates. Our
 findings highlight interactions among drivers of lichen decline (forestry,
 nitrogen, climate), functional traits (dispersal, lichen colour,
 sensitivity to nitrogen, water storage), and population processes
 (extinction/colonization).
description: 1. Description of methods used for collection/generation of
 data:  We extracted data from the Swedish National
 Forest Inventory (NFI) from two full inventory periods (IP; 1993-2002,
 IP1; 2003-2012, IP2) from plots in productive forest land. Formally
 protected forests were excluded. The design of the NFI includes
 stratification into five regions, with different sampling intensities and
 clustering of sample plots in square-formed tracts (4–8 plots per tract).
 The plots are circular with a radius of 10 m (314 m2) and are located
 around the tract perimeter (the length of tract side varies from 300 to
 1200 m among regions). Hair lichens are surveyed in permanent plots that
 are measured with a 10-year interval. Occurrence and
 maximum thallus length of <em>Alectoria</em>,
 <em>Bryoria</em> spp. and <em>Usnea </em>spp. were
 recorded on branches of <em>Picea abies </em>and the stem up
 to 5 m above ground. We selected eight variables of ecological importance
 for the studied lichens. Four were measured in NFI-plots: diameter at
 breast height (DBH) and crown limit (CRL) for each tree, whereas basal
 area (BAS) and stand age (AGE) were measured at plot level. We also
 extracted a landscape context variable (MAT, mature forest), defined as
 the percent of forest ≥60 years old, in a 100 m buffer around the centre
 of the plot. By using ArcGis version 10.3, we extracted MAT from SLU
 Forest map (SLU, 2018), where forest age was estimated in 25 m × 25 m
 cells by combining satellite and NFI data (Reese et al., 2003). MAT
 represents only one year (2005) as there were gaps in older
 data. Gridded climate data (4 km × 4 km) were obtained
 from the Swedish Meteorological and Hydrological Institute (SMHI, 2006;
 https://www.smhi.se/en/climate). We calculated mean annual temperature
 (TEMP) and mean total rain per year (RAIN) for each IP in all NFI-plots.
 RAIN was defined as the sum of precipitation in days with mean temperature
 ≥0°C, during which lichens can be active. We extracted deposition of
 atmospheric inorganic N (dry plus wet deposition; NDEP) from gridded data
 (20 km × 20 km) based on the Match model (Robertsson et al., 1999;
 www.smhi.se/klimatdata/miljo/atmosfarskemi). Mean annual N deposition was
 calculated for 1998−2002 and 2003−2012. 2. Methods for
 processing the data:  Lichen occurrence and thallus
 length We first estimated the total number of live
 <em>Picea </em>≥150 mm in each region and IP under a two-phase
 sampling design by taking stratum area, clustering in tracts, area of
 sample plots, and design weights for the sample trees into account. We
 then calculated the occurrence proportion of each lichen on the sample
 trees as the ratio between the estimated number of trees with presence of
 the lichen and the estimated total number of trees. The occurrence
 proportion was estimated for all region and IP combinations, across all
 regions, and separately for remeasured and new trees in IP2, together with
 corresponding 95% confidence intervals (CI). The estimates are
 approximately unbiased due to the large sample size. The mean thallus
 length of each lichen on occupied sample trees and 95% CI were estimated
 based on the same principles as described above. 
 Explanatory variables in all plots We
 calculated summary statistics for the explanatory variables across all
 NFI-plots with sample trees in each IP to evaluate the magnitude and
 direction of changes over time in these variables. Means and 95% CIs were
 calculated for all region and IP combinations, and across all
 regions.  Explanatory variables in plots with
 remeasured trees The trees that were remeasured (~50%)
 allowed us to examine how changes in lichen occurrence correlate with
 changes in variables over time in the plots. A lichen is either present
 (P) or absent (A) on a sample tree in each IP, and thus there are four
 occurrence trajectories (outcomes): persistence (PP), absence (AA),
 colonization (AP, ‘gain’) and extinction (PA, ‘loss’). We used Chi-square
 tests to examine the association between trajectories and regions for each
 lichen. Extinction and colonization rates were calculated following Yalcin
 & Leroux (2018). The extinction rate was calculated as the ratio
 of the number of trees where the lichen had disappeared in IP2 divided by
 the number of trees where it was present in IP1. The colonization rate
 represented the ratio of the number of colonized trees in IP2 divided by
 the number of trees where the lichen was absent in IP1. We also
 constructed maps to depict patterns of extinction and
 colonization..  We calculated the correlation
 coefficient (<em>r</em>) between all variable pairs to
 identify potential associations among the variables. We then used
 multinomial logistic regression (Hosmer et al., 2013) to identify the
 variables that best correlated with the occurrence trajectories and to
 examine the relationships. In such ‘multi-outcome’ models the dependent
 variable has more than two levels, in our case four levels. We fitted
 ‘single variable’ logistic regression models for each of the seven
 variables measured in IP1, the corresponding difference variables, and
 MAT, as well as ‘multiple variable’ models. A single variable model
 includes one explanatory variable, but it may appear in several
 transformations. For model-building we used a Multivariable Fractional
 Polynomials (MFP) approach (Hosmer et al., 2013) which finds non-linear
 transformations if sufficiently supported by the data, and removes weakly
 influential covariates by backward elimination (Sauerbrei &
 Royston, 1999). For this purpose, we used the R package mfp (Ambler
 & Benner, 2015). Although this package does not support
 multinomial logistic models, the suggestion by Hosmer et al. (2013) is
 that an individualized fitting approach, where separate binary logistic
 models are fitted, is useful for finding suitable non-linear
 transformations. Without taking the sampling design into account, we used
 this approach for finding preliminary main-effects models. After this
 stage, the preliminary main-effect models were refitted taking the complex
 sampling design of the NFI into account (Supporting Information Appendix
 2; cf. Ekström et al., 2018), after which we considered possible
 interactions between the main effects. Any model selected was considered
 preliminary until we evaluated its fit. We used the Akaike information
 criterion for the final model selection. Model performance was evaluated
 by calculating McFadden’s pseudo <em>R<sup>2</sup>
 </em>(Menard, 2000). The multiple variable multinomial logistic
 models were interpreted with odds ratios. 
description: For variable AGE1, stand age weighed by basal area in IP1, year, 
 -2 is a missing value for one sample tree