Weak recovery of epiphytic lichen communities in Sweden over 20 years of rapid air pollution decline

Abstract Epiphytic lichens are sensitive to deteriorating air quality, but levels of nitrogen and especially sulphur deposition have been in decline over most of Europe in recent decades. We assessed the response of epiphytic lichens to this decline, using data from long-term monitoring sites in Sweden. We analyzed 20 years of data to investigate temporal trends in lichen communities’ sensitivity to sulphur, nitrogen preference, species richness and alpha and beta diversity. We found only limited and partial evidence of recovery in the area that previously had high levels of deposition, and a decline in mean sulphur sensitivity at a northern site with low deposition levels throughout the monitoring period. The slow recolonization of sensitive species, even where environmental conditions are now suitable, is probably a result of impoverished regional species pools and the inherent limited dispersal capacity of many lichen species. We suggest due consideration of these factors in the use of epiphytic lichens as environmental indicators in a period of improving air quality.


Introduction
Lichens have long been recognized as good indicators of air quality, with perhaps the first attempt to quantify the relationship between lichen flora and pollution dating back to 1912 (Sernander 1912). Many lichen species are sensitive to air pollution and react rapidly, especially to sulphur and nitrogen, while some are more tolerant (Skye 1979). A thallus surface without protection, in combination with a non-specific uptake of mineral nutrients, makes lichens inherently vulnerable to pollution effects. This sensitivity is augmented by a slow growth rate, growth on substrata often exposed to air pollution, and an ability to absorb more sulphur dioxide (SO 2 ) at a given concentration than most vascular plants (Nash & Gries 1991). As a consequence, many lichen species have been adversely affected by the widespread high levels of sulphur deposition originating from the burning of fossil fuels in industrialized areas, which also spread further afield via long-range atmospheric transport (Gilbert 1986;Richardson 1988).
While this acid rain was for a long period a serious environmental concern, expanding controls on sulphur emissions (S) in many countries during the 1970s and 1980s, mainly through the UN Convention on Long Range Transboundary Air Pollution (CLRTAP; https://www.unece.org/env/lrtap/welcome. html.html), resulted in steep drops in deposition levels across most of Europe, an achievement that is considered an example of effective environmental legislation and international cooperation (Grennfelt et al. 2020). One outcome of these controls has been improved conditions for many lichen species. Although there are numerous studies showing the suitability of lichens as indicators of declining air quality, there are relatively few studies investigating the recovery of the lichen community as an effect of improved quality. In UK-wide data from 1970-2015 for example, an analysis of occupancy trends showed an increase in the combined bryophyte/lichen grouping of 36%, a change attributed to improvements in air quality (Outhwaite et al. 2020). Long-term monitoring data have shown a recovery of SO 2 sensitive species in areas where SO 2 deposition levels have trended strongly downwards in recent decades (Pescott et al. 2015), but failures in recolonization have also been recorded in some species while others increased (Bates et al. 2001). Earlier studies in London and surrounding areas show a similarly mixed response: some sensitive species extinct in 1970 had recolonized areas of London by 1980 (Rose & Hawksworth 1981) in what the authors characterized as a slight improvement, while few significant changes in the lichen flora were found on oaks in previously highly polluted areas from 1979-1990(Bates et al. 1990).
Attempts to reduce levels of nitrogen deposition (principally nitrogen oxides NO x from the combustion of fossil fuels and ammonium NH 4 from agriculture, hereafter referred to as N) have met with more limited success (Sutton et al. 2011). One result of this is that changes seen in recent decades in air pollution patterns do not always have a straightforward effect on lichens, since the much steeper decline generally seen in deposition of SO 2 than in N must be taken into account. Areas previously heavily affected by acidification can show a recovery of acid-sensitive species of epiphytic lichens, but simultaneously demonstrate a eutrophication effect due to persistently high deposition of nitrogen while the sulphur deposition has decreased (Friedel & Müller 2004). Long distance N deposition can adversely affect acidophytic species (van Herk et al. 2003), meaning that declines in acidophytic species may not always be a consequence of declines in S deposition alone. Experimental N addition suggests that over several years low levels of N deposition cause sensitive acidophytic species to be replaced by nitrophytic species (Munzi et al. 2010). One attempt to quantify the relative contributions of S and N to changes in lichen communities showed that combined pollutant deposition explained a majority of the variation in macrolichen abundance seen, with the respective contributions to the variation of 34.3% N and 22.4% S, indicating that at least in some areas N now has a greater impact (Giordani et al. 2014).
Cover of epiphytic lichens on mature living tree trunks is one variable in the ecosystem effect monitoring under the CLRTAP, that has been monitored every fifth year since 1997 at four permanent monitoring sites in Sweden (Fig. 1), along with monthly measurements of deposition levels of both nitrogen and sulphur compounds. Together with the strong gradient in deposition levels between the sites, from high levels on the south-west coast to very low levels in the north, this dataset provides a basis for investigating the long-term impacts of atmospheric deposition of pollutants on lichen communities. The monitoring sites ( Fig. 1) in this on-going study are located in long-term protected spruce forest nature reserves. Avoiding the confounding effects of forestry is particularly important for lichen studies as forest continuity is important for many species, due to the higher substratum quality found there and sufficient time available for slow colonization processes to occur (Fritz et al. 2008).
In this study, we use a 20-year time series (starting in 1997) on epiphytic lichens collected along a depositional gradient with decreasing deposition over time, to assess the extent to which the lichen community has recovered as a result of decreasing levels of air pollution. We hypothesize that lichen communities in the most polluted sites were depleted at the start of the monitoring period and will show a recovery during the studied years, while the lichen community in the pristine northern site should not show any temporal trends. However, as levels of S deposition have declined more sharply than N deposition, we expect the mean S sensitivity of the lichen community to increase at the polluted sites while the mean nitrogen preference at those sites will show more limited or no decrease.

Methods
The monitoring sites ( Fig. 1) vary in size from 3.7 to 45 ha. At each site, circular monitoring plots of 314 m 2 are located at the intersections of a 50 × 50 m grid (100 × 100 m at the largest site). Lichen monitoring is conducted on five randomly selected mature trees (diameter breast height 40-80 cm) from four randomly selected monitoring plots. However, due to severe storm and bark beetle disturbance in 2005-2008 at the Aneboda site (Weldon & Grandin 2019) the trees selected for lichen monitoring are restricted to those monitoring plots that still have living trees. The lichen monitoring follows a repeated measure design with inventories of the same trees every fifth year, resulting in four inventories at each site, as of 2020.
On each selected trunk, all epiphytic lichens are recorded using the point frequency method. A 40 × 40 cm transparent film with 400 regularly spaced holes is attached to the trunk at 120 cm above the ground, and the number of holes covering each lichen species under the film are counted.
Chemical concentrations at the sites are taken from data on wet throughfall deposition, which is most relevant for the availability of N/S to organisms under a forest canopy, and SO 4 is reported rather than SO 2 . However, data from Sweden-wide monitoring (Ferm et al. 2019) clearly demonstrate that SO x and SO 4 have followed a very similar trend of steep decline over the period in question. Funnel collectors are placed in sufficient numbers to permit measurements with at least 20% error of the mean with a 90% confidence interval for the measured ions, and samples are collected and processed at least monthly. Full details of the throughfall deposition sampling methodology can be found in part XIV of the ICP Forests manual (ICP Forests 2016) and section 7.5 of the Manual for Integrated Monitoring (1998). More information on all aspects of the monitoring programme can also be found in these documents. Further details of the sites can be found at DEIMS-SDR (Dynamic Ecological Information Management System -Site and dataset registry; https://deims.org/) using the permanent site identifiers given in the Supplementary Material (available online).

Statistical methods
For each tree, a lichen community air pollution sensitivity index was calculated as a weighted community mean value of the air pollution sensitivity values provided by Hultengren et al. (1991). The Hultengren sensitivity values range from 0 to 9 and are mainly an indicator of sensitivity to SO 2 . The higher the value, the less tolerant to acidity (i.e. the better air quality). Similarly, species preference for nitrogen was calculated based on the values given in Wirth (2010), with values ranging from oligotrophic (1) to eutrophic (9). Shannon diversity index values for each site/year combination were calculated using the R package vegan (Oksanen et al. 2019). Beta diversity between the communities on sample trees was investigated within each site using betapart in R (Baselga et al. 2018). Changes in beta diversity are decomposed into a turnover and a nestedness component, quantifying the extent to which changes are driven by species replacement or community homogenization (Baselga 2010). The analysis is based on Sørensen's dissimilarity index, using presence/absence data. We used the 'beta.multi' function to obtain a total value for turnover, nestedness and the overall beta diversity, measured as Sørensen dissimilarity for each site/year combination. We then used the 'beta.sample' function to resample these multiple plot dissimilarities 100 times and calculate the probability of a higher/lower value of each beta diversity measure in survey 4 relative to survey 1, yielding a P-value (Baselga 2010;Baselga et al. 2015). To ensure robustness of the results we also applied a complementary analysis of the pairwise matrices generated by the 'beta.pair' function, using a Mantel test (in vegan) to test for similarity between the results for survey 1 and survey 4. Note that as this is a test for similarity, we would generally expect results opposite to the resampling test for differences.
Changes in S-and N-sensitive species were analyzed, using the following definitions: N-sensitive species are those in the 1st quartile in the distribution of N preference values of all species found across all sites; S-sensitive species are those in the 4th quartile in the distribution of the S sensitivity values of all species found across all sites.
Temporal trends in diversity and the lichen community indices were assessed using a mixed model ('lme' in R package nlme (Pinheiro et al. 2019)) to account for the nested nature of observations (i.e. tree nested within plot as a random factor) and a first order autocorrelation structure to compensate for repeated observations by assigning time as a continuous covariate ('corCAR1' in nlme). Trends in sulphur and nitrogen deposition were tested for significance using simple linear regression with the 'lm' function in R. All data processing and analyses were performed in R version 3.6.2 (R Core Team 2019).

Results
The randomly selected trees on which epiphytic lichens were recorded were of several species but Picea abies and Pinus sylvestris L. were the most common at all four sites (Table 1). Picea abies is the most common species at Aneboda and Kindla, while Pinus sylvestris is the most common at Gårdsjön and Gammtratten. Deciduous species are also present at all sites, at lower frequencies.
Deposition Gårdsjön, while still having the highest deposition levels of N and S, has seen steep drops in the concentrations of airborne pollutants found in throughfall monitoring over the 20 years included in the dataset. Levels of SO 4 , NH 4 and NO 3 all declined significantly (P < 0.05, Fig. 2). Kindla and Aneboda have also shown significant reductions in S deposition. NH 4 and NO 3 concentrations declined significantly at Kindla, but at Aneboda NO 3 levels remained flat and there was an increase in NH 4 . Gammtratten has had low levels of both N and S deposition throughout the monitoring period, but even here levels of SO 4 and NO 3 have declined significantly.

Sensitivity index
The Hultengren sensitivity index increased at Gårdsjön (P = 0.009) and decreased at Kindla (P = 0.03) and Gammtratten (P = 0.0004) (Fig. 3, Table 2). At Gammtratten the decline in sensitivity takes the form of an initial drop followed by no clear trend. Although the change at Aneboda is not significant, there is a clear change in the distribution of results following the disturbances there, beginning in 2005 (Fig. 3).

Richness of S-and N-sensitive species only
The richness of S-sensitive species (defined as the 1st quartile in the distribution of the Hultengren sensitivity values of all species found across all sites) decreased at Gammtratten (P = 0.028) and Aneboda (P = 0.0078) and showed no significant change at Kindla and Gårdsjön (Table 6). The richness of N-sensitive species (defined as the 4th quartile in the distribution of the N preference values of all species found across all sites) decreased at Aneboda (P = 0.0041) and showed no significant change elsewhere (Table 7).
The tested response variables showed mainly decreases at the southern sites of Gårdsjön and Aneboda, while the two less affected sites Kindla and Gammtratten showed few significant changes over time (Table 8).

Beta diversity
Further insight into community shifts can be gained from analyzing the changes over time in beta diversity, taking the communities present on individual trees as the unit of comparison (Table 9).
At Gammtratten, the northern 'pristine' site, all measures of beta diversity were stable, with no significant changes found between the first and last surveys. Kindla was also largely stable, but with a small decline in overall beta diversity on the resampling test (although also significant similarity on the Mantel test). At Aneboda, however, there is an increase in turnover, a decrease in nestedness and an increase in overall beta diversity. At Gårdsjön, the most polluted site, there was no overall change in beta diversity but this masked a decrease in turnover and an increase in nestedness.

Discussion
We hypothesized that the lichen community in the more polluted area would show a recovery during the period of the monitoring programme, while the lichen community in the low pollution northern area would not change significantly, and that mean S sensitivity would increase at the polluted sites while the mean nitrogen preference would show more limited or no decreases. The most polluted site (Gårdsjön) did show an increase in mean S sensitivity and a decline in mean N preference but also declines in species richness and Shannon diversity. The two sites with intermediate pollution levels demonstrated no improvements and declines on some measures, while even the supposedly pristine far northern site (Gammtratten) showed a decline in S sensitivity.
The Gårdsjön site, while still showing the highest deposition levels, has seen the largest drops in both S and N deposition. Despite this, the recovery of the lichen flora is weak. An increase in the mean index of sensitivity to S was found but no increase in the most sensitive species, whose richness remained very low throughout the monitoring period. Similarly, a decrease in mean preference for N was seen, but no increase in the richness of the most N-sensitive species. Analysis of beta diversity shows no significant overall change but this masks a decrease in turnover and an increase in nestedness in the most recent survey, which could indicate increased homogenization of the community. Turnover, however, was still the dominant process in beta diversity changes, indicating a role for recolonization and/or local extinction although this clearly involves the more nitrophytic and acid-tolerant species rather than a return of sensitive species. Overall species richness and Shannon diversity index values have continued to decline over the monitoring period, which suggests that recovery of the lichen community is still far from complete after 20 years.
Interpretation of the Aneboda results is complicated by the extensive disturbances caused by storm and bark beetle damage at the site (Weldon & Grandin 2019). While there is an apparent large decline in the Hultengren sensitivity index and increase in N preference, these were not significant (although the large increase in variation may well be the reason for this lack of significance). Richness and Shannon diversity values, however, declined after the disturbances which began in 2005, and a decrease in nestedness and increase in turnover and overall beta diversity can be seen. However, given that the disturbances resulted in extreme changes in light, temperature, moisture and nutrient regimes at the site, as well as a patchwork of newly opened areas and surviving closed canopy areas, it is difficult to identify the specific contribution of atmospheric pollutants to these changes.  The Kindla site in contrast has lower levels of deposition than Gårdsjön but the mean Hultengren sensitivity index has decreased, indicating continuing effects from S deposition despite relatively low deposition levels. Mean N value has increased, indicating a eutrophication effect from N deposition, again with relatively low levels of deposition. Despite these shifts in community composition, however, no significant changes in richness or Shannon diversity index values were found, while the beta diversity results were inconclusive.
The northern 'pristine' site, Gammtratten, has as expected the lowest levels of both S and N deposition throughout the    monitoring period but despite this, the Hultengren sensitivity index declined from a mean of 3.96 in 2000 to 3.09 in 2015, a change also seen in the decline in richness of the most sensitive species. It is apparent (Fig. 3) that this can be characterized as a failure to recover from losses early in the monitoring period. No change in N preference at Gammtratten was found, indicating an absence of eutrophication effects, and no clear pattern of change can be seen in richness, alpha diversity or the beta diversity components, supporting our assumption that Gammtratten is to a large extent a pristine site with low previous and present deposition and hence minimal dynamics in the epiphytic lichen community. Nevertheless, the changes in S-sensitive species are concerning. Shannon Diversity −* −** n.s. n.s. * = P < 0.05, * * = P < 0.01. Minus sign indicates a significant decrease, while a plus sign indicates a significant increase. 'n.s.' indicates a non-significant change. Table 9. Total beta diversity, split into nestedness and turnover, for epiphytic lichens at monitoring sites in Sweden. The significance tests compare changes between the first and last survey. Sites are ordered from north to south. Beta diversity between the communities on sample trees was investigated with betapart in R (Baselga et al. 2018), where 'beta.SIM' was the measure of species turnover, 'beta.SNE' of nestedness and 'beta.SOR' was the total beta diversity, produced in all cases by the 'beta.multi' function. The resample test is based on using 'beta.sample' to resample the multiple location dissimilarities 100 times and calculate the probability of finding a higher/lower value in Survey 4 than in Survey 1. The Mantel test is applied to pairwise dissimilarity matrices generated by the 'beta.pair' function for each survey and is a test of similarity (i.e. we would expect opposite results to the resample test).
Lichen communities can respond quite rapidly to changes in air pollution levels, and inter-year changes can be detected (Loppi et al. 2004). In the 20-year period of this study we can indeed see changes in the lichen communities at the studied sites but the hypothesized recovery of sensitive species with improvements in pollution levels was only partially evident. At the most polluted site, steep declines in deposition of S and N have allowed an improvement in mean Hultengren sensitivity index and mean N preference. However, the least polluted site surprisingly showed a decline in S-sensitive species between the 2000 and 2005 surveys, from which it has not recovered. The species most sensitive to both S and N deposition can be lost at relatively low levels of deposition and cumulative deposition can be more important than current levels, meaning long-term perspectives are needed to understand the changes observed (Cleavitt et al. 2015). While both nitrogen and sulphur deposition have declined strongly, levels of both remain above estimates of preindustrial levels and may still be sufficient to negatively impact the most sensitive species. Modelling suggests that sulphur is still 30% above 1900 levels (which were almost certainly already elevated) and at current rates the decrease from the 1980s peak will continue until at least 2050, while reduced nitrogen deposition is around double that of 1900 and is unlikely to return to pre-industrial levels in the foreseeable future .
Lichen communities are also known to be very slow to fully recover from some disturbances. Studies analyzing community composition after fire, for example, have found that sites with the longest time since the last fire (several hundred years) show the greatest richness (Boudreault et al. 2000), while communities that are sensitive to disturbance can require up to 300 years to fully recover (Halpern & Spies 1995). While fire can be an extreme disturbance, fires very rarely have an impact over the same spatial or temporal scales as the long-distance transport of atmospheric pollutants. Disturbances over wide geographical scales and long periods can result in a regional species pool which is depleted of sensitive species (Cornell & Harrison 2014). A disturbance such as acid rain clearly does not have an influence only on a monitoring site but also on a much wider surrounding area. It is likely that some sensitive species were lost before the monitoring programme began, given that European emissions of sulphur dioxide peaked c. 1980 and emissions of nitrogen oxides c. 1990 (Grennfelt et al. 2020). While improved conditions imply that sensitive species should recolonize, this requires those sensitive species to be still sufficiently abundant in the regional species pool after this long period of elevated deposition of atmospheric pollutants (i.e. that there are local source areas from which to recolonize). If the regional species pool is depleted, recolonization must occur across much longer distances than a strictly local disturbance would involve, and the colonization of new sites is often a slow process in lichens (Dettki et al. 2000;Sillett et al. 2000;Buckley 2011). An example in Scandinavia is a Norwegian study of epiphytic lichen dispersal across isolated stands in a 170 km 2 area which found the first colonizing species took 40-50 years to appear and 100-150 years for the full community to establish (Gjerde et al. 2012). Another factor potentially limiting recolonization is that surrounding forested areas are mostly actively managed for timber production (as is the majority of all forest in Sweden). This results in a much larger proportion of young forest than would naturally be found, which restricts the colonization of those species which are more commonly (or even exclusively) found in old-growth forests (Lie et al. 2009). The combined effect of air pollution and forestry practices can result in the survival of isolated populations of sensitive lichens which are vulnerable to local extinction, with the nearest population source too far away for recolonization to occur. While long-term systematic monitoring of complete lichen communities is performed only at a small number of sites, by focusing on a single species we can see these trends at a larger scale. For example, a resurvey of 12 forest sites in southern Sweden after nine years found the sensitive macrolichen Lobaria pulmonaria (L.) Hoffm. present only at those sites where it was found in the first survey, despite suitable habitat being available within a few kilometres (a mean dispersal distance of only 35 m was achieved between surveys) (Öckinger et al. 2005). Another study in the same area found that at 13% of 66 sites monitored, L. pulmonaria had become locally extinct within 10 years, despite improving air quality (Öckinger & Nilsson 2010). For further discussion of L. pulmonaria and the regional species pool, see Supplementary Material (available online). While our focus here is on lichen diversity and community recovery, previous research focused on lichens as bioindicators has demonstrated that recolonization is often more heterogenous and complex than the declines seen under high levels of pollution. As early as 1989 the phenomena of 'zone skippers' and 'zone dawdlers' were documented, species whose presence or absence under conditions of rapidly improving air quality does not correspond with pollution levels as expected, but depends to a large extent on their dispersal abilities (Hawksworth & McManus 1989;Gilbert 1992). These results clearly show the importance of a patchwork of forested nature reserves as refuges and stepping stones in the preservation and recolonization of the depleted lichen community in landscapes affected by large-scale air pollution.
While we expected to see a clear improvement in indicators of lichen community health over a period of declining depositions, our results show a recovery that is at best partial, alongside further declines in the Hultengren sensitivity index at two sites and some indications of eutrophication and community homogenization. There are two likely explanations for this, the first being continued/cumulative deposition adversely affecting some sensitive species. A second factor likely to be important in explaining the observed limited recovery where conditions have improved is dispersal limitation from a depleted regional species pool adversely affecting recolonization. With continued declines in emissions, an eventual recovery is to be expected. However, given the wide geographical and temporal impact of the disturbance and the dispersal limitation of epiphytic lichen species, a full recovery of the pre-disturbance lichen community could take many decades. Our results confirm that lichens are sensitive to air pollution. However, the use of lichens as indicators during recovery from air pollution might be less reliable in cases where the regional species pool has been depleted (whether by pollution, management practices or other factors). Lichens may be good indicators of recovery from point source pollutions, but not necessarily in cases of large-scale transboundary air pollution, where the potential impact of regional species pools should also be taken into consideration.