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 >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).

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 Alectoria, Bryoria spp. and Usnea spp. were recorded on branches of Picea abies 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 Picea ≥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 (r) 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 R² (Menard, 2000). The multiple variable multinomial logistic models were interpreted with odds ratios. 

For variable AGE1, stand age weighed by basal area in IP1, year,  -2 is a missing value for one sample tree