Multiple drivers of large‐scale lichen decline in boreal forest canopies

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


| INTRODUC TI ON
The circumboreal forest constitutes 27% of global forest cover (Hansen et al., 2010). It is characterised by short growing seasons and cold winters with several months of snow cover, to which organisms have evolved various adaptations (Boonstra et al., 2016).
Both natural (fire, wind, insects) and human-induced disturbances (forestry, other land use) shape the structure, dynamics and function of the boreal forest, with about two-thirds under some type of management, mainly for industrial wood production (Gauthier et al., 2015). These forests are also under strong pressure from environmental hazards, including habitat fragmentation, air pollution and climate change, with severe consequences for biodiversity and ecosystem services (Esseen et al., 1997;Moen et al., 2014;Pohjanmies et al., 2017). Global warming is stronger at the high latitudes of the boreal biome than in southern biomes, implying multiple risks to forest organisms (IPCC, 2013;Venäläinen et al., 2020).
Lichens are common throughout the boreal biome and form rather homogeneous communities on the forest floor and in canopies (Ahti, 1977;Payette & Delwaide, 2018). This study focuses on hair lichens within the filamentous genera Alectoria, Bryoria and Usnea, which are dominant epiphytic species in boreal forests . These lichens show a clear vertical zonation, with dark Bryoria in the upper canopy and pale Alectoria and Usnea in the lower canopy (Coxson & Coyle, 2003). This zonation is mainly shaped by their colour and the function of species-specific sunscreening fungal pigments in the cortex (Färber et al., 2014). Niche differentiation is also linked to anatomical and morphological traits that influence uptake, storage and loss of water (Esseen et al., , 2017Gauslaa, 2014). Hair lichens dominate canopies in unmanaged boreal forests (Boudreault et al., 2015;Dettki, & Esseen, 1998). In fact, these lichens were so abundant that a century ago they were considered a serious problem for forestry, erroneously assumed to reduce tree growth, or even kill trees (Romell, 1922). Yet, these lichens do not significantly impair tree growth but instead provide important ecosystem functions, such as food for ungulates, particularly reindeer/caribou and small mammals, nesting material for birds, as well as microhabitat and food for invertebrates (Asplund & Wardle, 2017;Pettersson et al., 1995). Hair lichen abundance is also an easily visible indicator of general diversity of epiphytic lichenized and non-lichenized fungi. These lichens are particularly sensitive to air pollution (Geiser et al., 2019(Geiser et al., , 2021 and forestry (Esseen et al., 1996(Esseen et al., , 2016Hauck, 2011;Lesica et al., 1991;Stevenson & Coxson, 2007). Overall, however, we lack information about the speed and scale of changes in hair lichen populations and how different global change drivers interact and influence species with different functional traits.
Here, we have analysed how global change drivers influence changes in occurrence and length of hair lichens over 10 years in the lower canopy of Norway spruce (Picea abies (L.) Karst.; henceforth referred to as Picea) in managed forests in Sweden. Our analysis is based on >6000 trees surveyed 1993-2012 in the National Forest Inventory (NFI). We focus on changes in habitat quality by monitoring lichens when the host is present. This allowed us to examine how global change drivers interact and affect extinction and colonization processes. We hypothesized the following: (1) that lichen occurrence declines due to industrial forestry, driven by short logging cycles and an unfavourable microclimate (i.e. low light) in dense forests (Esseen et al., 2016); (2) that Alectoria, being associated with old forests (Esseen et al., 1996), declines faster than Bryoria and Usnea; (3) that lichen occurrence and thallus length (the length of the vegetative body) increase and recover in southerncentral regions in response to reduced anthropogenic N deposition since 1980 (Engardt et al., 2017). Alectoria and Bryoria are more sensitive to N deposition than Usnea (McCune & Geiser, 2009) and should respond faster to N reductions; (4) that climate change has genus-specific effects on lichen occurrence. Here, a warmer and wetter climate, with more rain and less snow, should favour the occurrence of Alectoria and Usnea, as rainfall drives growth in these lichens (Phinney et al., 2021). In contrast, Bryoria, adapted to cold and dry climates with low rainfall (Esseen et al., 2017), should decrease in such climates.

| Study area
The study area covers the whole of Sweden and spans latitudes 55-69°N (c. 1500 km). Sweden has 279,000 km 2 of forest, of which 235,000 km 2 is productive (site productivity ≥1 m 3 ha −1 year −1 ; SLU, 2021). Ninety-four percent of the productive forest is managed while 5.7% is formally protected (SCB, 2020). A large proportion of the protected forest forms intact landscapes along the eastern slopes of the Scandinavian Mountain range .
The study area covers large gradients in climate, vegetation, and anthropogenic impact. The climate ranges from cold temperate moist with snow in the north to warm temperate moist in the south, with warm summers in most of the country. The boreal zone covers most of the study area ( Figure 1a) and is dominated by coniferous trees (Picea, Pinus sylvestris L.). The hemiboreal zone, a transition between the boreal and temperate zone, covers most of southern Sweden.
The temperate zone forms a narrow belt in the south and southwest.
The current land use in Sweden is dominated by forestry (58% of the land area), followed by agriculture (8%; SLU, 2021). Reindeer husbandry is practiced in ~50% of the country and overlaps with forestry, which results in conflict of interest between traditional indigenous land use and modern industrial forestry (Horstkotte & Moen, 2019). Commercial forestry on a larger scale started with some selective cutting of large trees in the south-central parts in the 18th century and harvest of trees for charcoal in the region of Bergslagen, central Sweden, followed by a northward expanding timber frontier during the 19th century (Östlund & Norstedt, 2021;Östlund et al., 1997). The high-grading of large timber trees in old-growth and multi-layered forests was followed by a period of uneven-aged forestry and was then gradually replaced by even-aged forestry at around the mid-20th century (Lundmark et al., 2013). Current forestry is dominated by clear-cutting followed by planting of conifer monocultures, which are regularly thinned. Rotation cycles range from c. 50 years in the south to c. 120 years in the north (Fries et al., 2015). Eighty-two percent of the productive forest is ≤100 years old (SLU, 2021). Picea dominates by volume (41%), followed by P. sylvestris (39%) and Betula spp. (12%). The growing stock has increased with 106% since the 1920s, at which time the growing stock had decreased during the late 19th century. This trend mainly results from the efficient, production-oriented forestry, including cutting of unproductive stands, ditching, thinning, N fertilization and planting of seedlings on clear-cuts (Östlund et al., 1997). Global warming, N deposition and elevated CO 2 concentration have probably also enhanced forest growth (Elfving & Tegnhammar, 1996;Venäläinen et al., 2020).
The southern and central regions have earlier been exposed to high levels of SO 2 pollution, which has substantially decreased since 1970 (Vestreng et al., 2007). N deposition has currently greater impact on forests and exceeds 10 kg N ha −1 year −1 in southern regions (Pihl Karlsson et al., 2011), while deposition is low in northern regions.

| The Swedish NFI
We extracted data from the Swedish NFI, which is a broad-scale monitoring program that provides information about wood resources, ecosystem services, biodiversity, and carbon sequestration (Fridman et al., 2014). The design includes stratification into five F I G U R E 1 (a) Map of Sweden with five regions and National Forest Inventory plots clustered in tracts, with ≤4 (region 5) or ≤8 plots per tract (region 1-4). Region 1, northern boreal; 2, middle-northern boreal; 3, southern-middle boreal; 4, mainly hemiboreal; 5, temperate. (b) Alectoria sarmentosa (Witch's hair), a pale hair lichen associated with old coniferous forests in the boreal biome regions, with different sampling intensities and clustering of sample plots in square-formed tracts (4-8 plots per tract; Figure 1a).
The plots are circular with a radius of 10 m (314 m 2 ) and are located around the tract perimeter (the length of tract side varies from 300 to 1200 m among regions). About 200 forest, vegetation and site variables are recorded on each plot. The NFI has a systematic program for quality assurance, including training, calibration, and control inventory (Fridman et al., 2014). The data in this study were collected by ~150 field personnel.

| Lichen data
Hair lichens are surveyed by the NFI in permanent plots that are measured with a 10-year interval. We used NFI data from two full inventory periods (IP;1993-2002-2012 from plots in productive forest land (see above). Formally protected forests were excluded as they were not measured in IP1. Hence, our sample consists of managed forests with active forestry, but includes a small proportion of voluntarily set aside areas (~7%; SCB, 2020). We studied lichens on live Picea with DBH ≥150 mm and included all NFI plots with at least one such tree present. The sample tree was the first Picea encountered along the measurement direction in the plot. Therefore, a new tree could be selected in IP2, due to ingrowth, even if the old tree remained. The year and cause of death were recorded for trees that died or were cut. In such cases, a new tree was selected in IP2, while plots without sample trees were excluded. New plots were also included in IP2 when at least one Picea reached a DBH ≥150 mm. The sample included a total of 6140 trees (Table S1).
Occurrence and maximum thallus length (a proxy of lichen mass; Esseen, 2006;McCune, 1990) of Alectoria, Bryoria spp. and Usnea spp. were recorded on branches and the stem up to 5 m above ground. Our study is restricted to genus level as species are difficult to identify. However, Alectoria is only represented by A. sarmentosa (Ach.) Ach. (Figure 1b), a red-listed (near-threatened) species associated with old coniferous forests (Esseen et al., 1996)

| Explanatory variables
We selected eight variables of ecological importance for the studied lichens (Table 1). 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.

| 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 (Toet et al., 2007; Appendix S1).
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. We evaluated changes in regions and in the whole country by calculating 95% CIs for the difference in occurrence and thallus length between IPs. The change is significant (p < .05) if the CI do not overlap zero.

| Explanatory variables in all plots
We calculated summary statistics for the explanatory variables (Table 1) 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. However, MAT was only summarized by region as data represented 1 year (see above).

| 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 and 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. Some overestimations of colonization and extinction rates are likely, as short thalli may have been overlooked.
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 'multioutcome' 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 modelbuilding, we used a multivariable fractional polynomials 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) (Menard, 2000).
The multiple-variable multinomial logistic models were interpreted with odds ratios. Consider a model with a single continuous explanatory variable x and a response variable Y. Then, if j is one of the occurrence trajectories, the probability that Y = j is assumed to depend on the value of x and is denoted by Pr (Y = j|x). Absence was used as our baseline trajectory and is labeled as Y = 0. The odds ratio of trajectory Y = j (e.g. extinction) versus Y = 0 for a one-unit increase in the explanatory variable is then given by where Pr (Y = j|x) ∕ Pr (Y = 0|x) is the odds at x that the trajectory is j, given that it is either j or 0. If there are other explanatory variables than x present in the model, these are kept fixed when computing OR j . The odds ratio OR j is a measure of how much more likely or unlikely (in terms of odds) it is for occurrence trajectory j to be present among those trees with a one-unit increment in an explanatory variable x as compared to those with no increment in this variable, while holding the other explanatory variables fixed. The odds ratio is significantly different from 1 when the 95% confidence band does not overlap 1. The analyses were done with r version 4.0.3 (R Core Team, 2020).

| Host trees
The number of sampled Picea decreased with 6% over time (Table   S1). Of the trees, 52.4% in IP1 were remeasured, 17.5% were cut, 1.8% died and 28.4% were only measured once. The proportion of cut trees increased from 8.8% in region 1 to 21.0% in region 4.
Accordingly, the proportion of new trees increased along the same gradient. One fifth (20.8%) of the plots in IP2 were in forests subject to thinning.

| Lichen occurrence
The occurrence proportion of all lichens decreased over time.
The distributions of Bryoria and Usnea were more stable.
Remeasured trees had higher lichen occurrence than new trees in IP2, particularly in Alectoria, 0.218 versus 0.130, respectively (Table S2). Extinction and colonization rates were high in several regions (Table 2; Figure 3), indicating substantial turnover, especially in Alectoria and Usnea. The Chi-square test of association between occurrence trajectory and region was significant (p < .05) for all lichens, showing that trajectories varied by region.

| Thallus length
The length of Alectoria in IP1 (mean 18.1 cm across regions) was twice as high as Usnea (8.6 cm), with Bryoria intermediate (13.6 cm).

Length of Bryoria increased slightly over time in regions 3-5 and
Usnea in regions 2-4, whereas Alectoria did not change ( Figure S2).
However, the change in length across the country was only significant for Usnea.

| Explanatory variables in all plots
The individual NFI plots spanned a large range in TEMP (−1.7 to 8.6°C) and RAIN (335-1290 mm). All variables showed clear latitudinal trends from north (region 1) to south (region 5; Figure 4).

| Explanatory variables in plots with remeasured trees
The correlations were relatively weak between most variable pairs, but significant due to the large sample size (Table S3) (Table S4). The difference between extinction and colonization probabilities in Figure 5 indicates the direction and magnitude of changes in lichen occurrence across the range of each variable.

TA B L E 2
Occurrence trajectories of studied lichens over 10-year on remeasured trees included in multinomial logistic regression models. Values in parentheses represent extinction and colonization rates (see Section 2 for explanation F I G U R E 3 Distribution of extinction (filled red circles) and colonization (filled blue circles) of the studied lichens on remeasured Picea over a 10-year period (1993-2002; 2003-2012). N = 2196 The multiple-variable model for Bryoria had higher explanatory power (R 2 = .337) than models for Alectoria (R 2 = .218) and Usnea (R 2 = .186; Table S5). The models included five or six variables, but only two of these were difference variables (TEMP∆ and DBH∆).

| DISCUSS ION
Habitat loss, fragmentation and degradation are main drivers behind loss of biodiversity in forests across the globe (Fischer & Lindenmayer, 2007;Haddad et al., 2015). Here, we show that hair lichens decline in managed forests across broad spatial scales even when an important host tree (P. abies) in the boreal biome is present. Our key finding of significant declines in Sweden over 10 years is alarming, as these lichens already experienced F I G U R E 6 Odds ratio functions (solid lines, log-scale) for occurrence trajectories of Alectoria over 10 years on remeasured Picea with 95% confidence bands (dashed lines) when adding one unit to either TEMP1, TEMP∆, AGE1, RAIN1, MAT or DBH∆ while holding other variables constant in multiple variable multinomial logistic regression models. Note that AGE1 interacts with TEMP∆ and is displayed on the x-axis, while TEMP∆ is displayed on the x-axis for AGE1. Thin grey horizontal lines indicate an odds ratio of 1.0. Horizontal coloured lines below the x-axis indicate the intervals where the odds ratios of the trajectories are significantly different from 1.0. DBH, diameter at breast height; MAT, mature forest large-scale declines also before the 1990s (Bruteig, 1993;Esseen et al., 2016;Kuusinen et al., 1990). The speed of the Alectoria decline in region 3 (~50% loss) is remarkable, as the distribution of this lichen had once been centred in this region of Fennoscandia (Ahlner, 1948). We base our interpretation of causal mechanisms on state and change in variables influencing lichen occurrence on all trees, and on logistic regression models for occurrence trajectories on remeasured trees, where each tree was monitored over time.

| Forestry and forest structure
Our first hypothesis, predicting that industrial forestry causes a decline of hair lichens is supported by the following points: (1) Declines were observed in all regions, and forestry is the only driver with significant impact across the country.
(2) The estimated annual loss rates of the lichens (0.5%-1.7%) are comparable to that ~1% of the forest area is logged each year, mainly by clearcutting of older forests. Moreover, the proportion of logged forests older than 120 years has increased since 2000 (Fries et al., 2015). Logging of old, multi-layered forests with long continuity of cover (Svensson et al., 2019) is harmful for these lichens (Dettki & Esseen, 1998;Esseen et al., 1996). (3) The new trees in IP2 had fewer lichens than remeasured trees, showing that the transformation from naturally to artificially regenerated forests contributed to the decline. (4) Short rotation cycles are incompatible with accumulation of high standing crop of these lichens (Dettki & Esseen, 2003). (5) Alectoria, with the strongest association to old forests, showed the steepest decline, supporting both our first and second hypothesis. This is not explained by growth responses as Alectoria grows faster than Bryoria and Usnea in lower canopy (Phinney et al., 2021). Instead, dispersal limitation and low propagule availability contribute to the scarcity of Alectoria in production forests. This lichen disperses with large thallus fragments, which are less suitable for long-distance dispersal (Dettki et al., 2000;Esseen et al., 1996), whereas Bryoria and Usnea disperse effectively by numerous small fragments and symbiotic propagules (Clerc, 2011;Esseen, 1985), which rapidly colonise tree saplings (Hyvärinen et al., 1999). (6) MAT was included in the multiple-variable models for all lichens, suggesting that landscape fragmentation also influenced occurrence trajectories. The lower MAT in southern regions implies greater proximity to edges and increased risk for adverse effects on microclimate as well as higher N deposition (Harper et al., 2015; F I G U R E 7 Odds ratio functions (solid lines, log-scale) for occurrence trajectories of Usnea over 10 years on remeasured Picea with 95% confidence bands (dashed lines) when adding one unit to either TEMP1, NDEP1, MAT, DBH1 or DBH∆ while holding other variables constant in multiple variable multinomial logistic regression models. See Figure 6 for additional explanations. DBH, diameter at breast height; MAT, mature forest Schmidt et al., 2017). It also reflects less adjacent forests functioning as source of dispersal propagules (Bartemucci et al., 2022;Dettki et al., 2000).
The low explanatory power of forest variables is partly due to the overriding effects of the large latitudinal gradients in climate and N deposition. However, it is well established that dense, productive conifer forests with high basal area provide poor growth conditions for lichens (Esseen et al., 2016;Hauck, 2011). Productive forests with a high CRL have few large, live branches in the lower canopy. Dead branches may be numerous but are smaller than live branches. Hence, the increase in CRL in all regions suggests that substrate limitation contributed to the decline. Moreover, lichen growth in lower canopy is substantially limited not only by low light in dense stands, such as spruce plantations, but also by less frequent hydration, due to rain interception in upper canopies. Bryoria, with dark cortical pigments, suffers more from low light conditions than Alectoria and Usnea, with pale pigments, as the former has a higher light compensation point for photosynthesis (Coxson & Coyle, 2003;Gauslaa et al., 2020). Alectoria and Usnea are adapted to the shaded lower canopy while Bryoria is adapted to the exposed upper canopy and open forests (Esseen et al., 2017;Färber et al., 2014).
Successional changes as forests become older and darker cause an F I G U R E 8 Odds ratio functions (solid lines, log-scale) for occurrence trajectories of Bryoria over 10 years on remeasured Picea with 95% confidence bands (dashed lines) when adding one unit to either TEMP1, TEMP∆, MAT, NDEP1, AGE1 or CRL1 while holding other variables constant in multiple variable multinomial logistic regression models. Note that MAT interacts with NDEP1, which is why NDEP1 is displayed on the x-axis for odds ratios of MAT, whereas odds ratios of NDEP1 are displayed as contour plots with MAT on the y-axis. For NDEP1, all odds ratios are significant for persistence and extinction, while odds ratios for colonization are significant where MAT is below 74%. See Figure 6 for additional explanations. CRL, crown limit; DBH, diameter at breast height; MAT, mature forest upward vertical displacement of Bryoria and may partly explain declines on remeasured trees.

| Nitrogen
The decline on remeasured trees and similarly recorded declines in protected forests (Jönsson et al., 2017) show that other drivers also affected the lichens. Although thallus length increased slightly in some regions, we found no increase in occurrence in southern-central regions in response to the reduced anthropogenic N deposition since ~1980 (Engardt et al., 2017). Moreover, the regression coefficients for NDEP∆ were not significant in the multiple variable models when other variables were included. Therefore, our third hypothesis predicting a recovery following the N reductions was not supported. Instead, we found strong evidence that high N deposition continues to limit occurrence of hair lichens. NDEP1 had significant explanatory power in the single variable models for all lichens as well as in the multiple variable models for Bryoria and Usnea. However, RAIN1 replaced NDEP1 in the multiple model for Alectoria. These variables were strongly correlated, and the odds ratio plots for RAIN1 in Alectoria and NDEP1 in Usnea were similar. The odds of extinction in Alectoria exceeded the odds of colonization at high rainfall (with high N), suggesting that N deposition contributed to its decline. High rainfall with a sufficient and not too high N concentration boosts growth of this lichen (Phinney et al., 2021).
High N deposition causes an imbalance between the photo-and mycobionts, loss of the integrity of the lichen symbiosis and increased risk for parasitic fungal attacks (Carter et al., 2017). The harmful effects of N also accumulate over time (Johansson et al., 2012), causing time lags in responses. Hair lichens are oligotrophic and N-sensitive and have critical loads (above which N is harmful) in the range of 2-6 kg N ha −1 year −1 (Esseen et al., 2016;Geiser et al., 2019Geiser et al., , 2021McCune & Geiser, 2009). Our data, based on the MATCH model, show that N deposition exceeded 5 kg N ha −1 year −1 in most of the south-central regions. However, N deposition in Picea forests can be up to 6 kg N ha −1 year −1 higher than estimates from this model (Karlsson et al., 2019). Hence, N deposition likely exceeded critical loads for all lichens in southern-central regions. Alectoria is more sensitive to N than Usnea (Esseen et al., 2016;McCune & Geiser, 2009) accelerating its decline.
N deposition contribute to the low occurrence of Bryoria in southerncentral regions but cannot explain its northern decline. Weldon and Grandin (2021) found only weak recovery of epiphytic lichens over the last 20 years in southern Sweden. Our data suggest >30 year recovery lag of hair lichens following the improved air quality, in accordance with studies from central Europe showing strong recovery delays in response to lower S and N deposition (Hauck et al., 2013;Schmitz et al., 2019).

| Climate
The occurrence trajectories correlated with TEMP1 in the models, reflecting the lichens macroclimate preferences along the latitudinal gradient (Esseen et al., 2016 The small increase in rain had marginal effect on lichen occurrence. However, precipitation is predicted to increase with 10%-40% during this century, particularly in northern regions (Sjökvist et al., 2015). Alectoria and Usnea grows better in wetter and warmer climates and may expand northwards unless light and dispersal are limiting. Alectoria grows faster than Usnea across continental to oceanic macroclimates (Phinney et al., 2021) due to a larger internal water storage pool (Eriksson et al., 2018;Esseen et al., 2015), extending hydration and carbon gain. Bryoria, in contrast, has large external water storage pool (Esseen et al., 2017) causing suprasaturation depression of photosynthesis, and in situ decomposition if exposed to prolonged hydration (Coxson & Coyle, 2003). This lichen may face periodic dieback in future wetter climates (Gauslaa, 2014;Goward, 1998), shrinking its distribution, unless forests become more open.
Overall, however, stochastic disturbances had less impact on lichens than anthropogenic stressors in the present study. Forest disturbances are predicted to increase in the boreal biome under future climates (Seidl et al., 2017;Venäläinen et al., 2020) and constitute potential threats to epiphytic lichens.

| Implications
The low occurrence and short length of the lichens, compared with natural forests (Esseen, 2019), implies that managed boreal forests have substantially lower standing crop than natural forests. This is reinforced by that we studied NFI plots with larger trees, which usually support higher lichen mass than smaller trees (Esseen et al., 1996). Clear-cuts and young, managed stands have none or low mass of these lichens (Dettki & Esseen, 1998 (Heggberget et al., 2002 Sandström et al., 2016). Loss of hair lichens will also affect other ungulates and small mammals and decrease abundance of canopy-living insects and spiders, constituting important prey for birds (Pettersson et al., 1995).
Sweden has had programmes for conservation of biodiversity in managed forests since the 1980s, focusing on retention (since the 1990s) and forest certification (Gustafsson & Perhans, 2010;Simonsson et al., 2015). Although retention of live and dead trees through the regeneration phase can mitigate some harmful effects of forestry on biodiversity (Fedrowitz et al., 2014;Rosenvald & Lõhmus, 2008), our results suggest that the current application of retention forestry in even-aged silviculture cannot support high standing crop of hair lichens. Instead, natural disturbance-based management strategies applied at multiple scales have large potential to maintain heterogeneity, microclimates, processes and biodiversity in forest ecosystems (Kuuluvainen et al., 2021;Messier et al., 2019). However, it is a challenge to integrate such multi-scale conservation approaches in the intensively managed production forest matrix (Felton et al., 2020). A spectrum of methods must be applied to secure substrates and suitable microclimates for hair lichens and other epiphytes under future climate scenarios (Ellis & Eaton, 2021), including: (1) uneven-aged, continuous-cover forestry, such as partial cutting (Stevenson & Coxson, 2007), (2) longer rotation cycles (Dettki & Esseen, 2003), (3) preservation of groups or large patches of lichen-rich trees to reduce edge influence (Esseen, 2019;Stevenson & Coxson, 2003), (4) variable density thinning to diversify substrates and microclimates (Kuuluvainen et al., 2021;Lehmkuhl, 2004) and (5) spatial planning to protect old, lichen-rich forests such that propagule availability is secured throughout landscapes (Dettki et al., 2000).

ACK N OWLED G EM ENTS
We thank the ~150 field personnel for collecting the NFI data and

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

AUTH O R CO NTR I B UTI O N S
P-AE designed the study and wrote the manuscript with help from all co-authors. BW extracted the NFI data and P-AE the other data, AG estimated lichen occurrence and length, ME performed the logistic regressions and P-AE performed the other analyses. ME, AG, GS and BW contributed to statistical aspects, while BGJ and KP contributed to ecological aspects. All authors contributed critically to the drafts and gave final approval for publication.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in Dryad at https://doi.org/10.5061/dryad.2ngf1 vhq5.