Index of Atmospheric Purity reflects the ecological conditions better than the environmental pollution in the Carpathian forests

The Index of Atmospheric Purity (IAP) is a popular tool used for the assessment of air quality in polluted urban areas, on the basis of phytosociological data of epiphytic lichen communities. We hypothesized that this indicator could also be used in less polluted forest areas to determine the quality of ecological conditions for lichens. The aim of the present study was to verify the use of IAP method in the assessment of environmental pollution, and alternatively for the assessment of general ecological conditions in protected mountain forests of Gorce National Park (Polish Western Carpathians) based on the epiphytic lichen biota associated with Picea abies. The spatial distribution of IAP values on monitored sites in GNP was compared with: 1) spatial distribution of accumulated sulfur, nitrogen, selected heavy metals, and total heavy metals in Hypogymnia physodes thalli in 1993 and 2018 (30 sites), 2) mean ecological indicator values characterizing species requirements for light (L), substrate reaction (R) and nutrients (N), in 1993, 2013 and 2018 (33 sites). Generalized linear model and redundancy analysis were performed for disclosing most influencing factors affecting lichen communities. The study revealed a few negligible relationships between IAP values and accumulation of such elements as Ni, Mn, Cd, and Cr in both monitoring periods. Simultaneously, IAP can be useful for the identification of forest areas with a high degree of naturalness.


Introduction
Due to water and toxic elements absorption directly from the atmosphere, lichens are regarded as excellent bioindicators of air pollution (Gries 2008;Sujetovienė 2015;Will-Wolf et al. 2015). Analytical techniques used for air quality assessment as well as regular monitoring of urban areas allow for the confrontation of the results from biological monitoring and the instrumental measurements of air pollution. Results of studies conducted in conurbations or near point sources of pollution usually confirm the legitimacy of using lichens in the bioindication of air pollution in these areas, for example those caused by sulfur dioxide (Białońska and Dayan 2005;Orlova et al. 2015;Mateos and González 2016), nitrogen oxides (Gaio-Oliveira et al. 2005;Tretiach et al. 2007;Mateos and González 2016) and heavy metals (Kularatne and De Freitas 2013;Parzych et al. 2016). The accumulation of toxic elements in lichen thalli is negatively correlated with an increasing distance from the main sources of pollution (Salemaa et al. 2004;Attanayaka and Wijeyaratne 2013). This correlation is clearly seen in the urban landscape, as air pollution is usually the main ecological factor determining the composition of lichen biota (Seaward 2015). For areas characterised by higher heterogeneity or degree of naturalness, the interpretation of an effect of urban and industrial pollutants on the environment is more complicated (Poikolainen et al. 2000;Gombert et al. 2004;Giordani 2007;Agnan et al. 2017). This is caused by the synergistic or antagonistic effect of many ecological factors (e.g. climate, topography, habitat, anthropogenic impact) determining the sensitivity of individual lichen species to the same concentrations of harmful substances (e.g. Van Dobben et al. 2001;Brunialti et al. 2010;Liu et al. 2016).
The Index of Atmospheric Purity (IAP) relies on phytosociological indicators and is a popular tool used for the assessment of air quality in urban areas. This method assumes that the growth of epiphytic lichens in urban areas is primarily determined by air quality. The basic formula developed by LeBlanc and De Sloover (1970) has been modified throughout the years. For example, Kricke and Loppi (2002) compared twenty distributions based on various versions of the IAP with the results of instrumental measurements conducted in the late 1980's (Ammann et al. 1987). The highest coefficient of determination (R 2 = 0.97, p < 0.05) was obtained by applying the formula based on the sum of frequencies of species found on a site. This variant of IAP was later used by Gombert et al. (2004) while conducting research in the urban area of Grenoble, France. Researchers concluded that the obtained values of IAP were more accurate when considered as a degree of anthropogenic transformation of the environment than as a mean annual concentration of air pollutants. Sites with high IAP values were dominated mainly by species without any strong preference for specific pH of the bark and trophic conditions, while nitrophilous and acidophilous species prevailed on sites characterised by low IAP values. Gombert et al. (2004) found no relationships between the distributions of the modified IAP values and mean annual concentrations of both sulfur and nitrogen oxides in the air as well as the distance from the main sources of pollution. This means that the use of the bioindication method based on the IAP criteria has to be adapted to local conditions and lichen biota and requires regular validation regardless of the degree of urbanisation. Despite the IAP is designed to describe specifically the urban areas, some authors have put effort into explaining the IAP values obtained from air-polluted forests (e.g. Poikolainen et al. 2000;Jeran et al. 2002;Jayalal et al. 2017), which can be considered wrong.
The doubts surrounding the IAP method application in air quality assessment does not preclude its use in bioindication of the forest environment, i.e. the quality assessment of ecological conditions for lichens. Conclusions resulting from this kind of approach would have to be general, as in more heterogeneous ecosystems the growth of lichens is affected by a more diverse range of limiting ecological factors, such as: insolation, relative humidity, competition for space, the age and species of trees, etc. (e.g. Gibson et al. 2013). With regard to mountain forests, an inference is even more difficult, as demonstrated studies carried out in Carpathians by Bytnerowicz et al. (2002) have shown that the relationship between pollution and environmental factors can change on a daily, seasonal and multiannual scale. For this reason, a particularly long-term research on the application of the IAP method in mountain forests is recommended to identify these environmental factors that determine the dynamics of lichen communities.
The aim of the present study was to evaluate the suitability of the IAP, as originally proposed by LeBlanc and De Sloover (1970), as follows: 1) as a bioindicator-based assessment of pollution in large areas of natural forests in the Western Carpathians caused by toxic elements and to compare the obtained results with the concentrations of these elements instrumentally measured in Hypogymnia physodes thalli; 2) as an alternative use in an assessment of general ecological conditions in the dynamically changing forest ecosystems of the Carpathians. The following working hypotheses were adopted: 1) there is a consistent lack of significant (α = 0.05) relationships between IAP values and the concentrations of elements measured separately in two study years -1993 and 2018; 2) IAP values for study areas correspond with ecological indicator values (Wirth 2010) in each of the three study years: 1993, 2013 and 2018.

Study area
The study was conducted in Gorce National .850" E) located in the Gorce range in the Polish part of the Western Carpathians. The Gorce ranges are flysch, mediumheight mountains with the highest peak, Turbacz, at 1311 m a.s.l. (Cieszkowski 2006).
Gorce National Park was established mainly to protect remnants of old-growth Carpathian forest with preserved patches representing different variants of the lower belt beech forest Dentario glandulosae-Fagetum and fir-spruce forest Abieti-Piceetum montanum as well as subalpine spruce forest Plagiothecio-Piceetum (Loch and Armatys 2008). The subalpine zone ranges above ca. 1050 m a.s.l. Fragments of primeval forest cover only 1.26% of the total forested area of GNP, which is 6313.6 ha. Natural forests and managed forests of a semi-natural character (with species composition consistent with the habitat, mainly formed by natural regeneration) account for 90% of this area (as of 1997; Loch and Armatys 2008).
The climate of Gorce range is mainly influenced by montane conditions. The diversified topography defines distinct differences in temperature, insolation, humidity, wind speed and the persistence of snow cover between sites that are not located far apart as well as for significant local fluctuations in these parameters. The mean annual temperature in the foothills is 6°C and drops to 3°C on the top of Turbacz Mt. The total annual precipitation ranges between 800-1200 mm and varies depending on the altitude above sea level. Most of the pollutants are carried together with the wind, usually blowing from the N and NW directions, and originate in the distant industrial region of Upper Silesia and Kraków located about 60 to 100 km northwest and north of the Gorce range. The wind direction is also locally determined by the course of river valleys. There are also periodic foehn winds that reach high speeds causing numerous windthrows in large areas, accelerating decomposition of spruce stands (Miczyński 2015).

Field data collection
The study was conducted in three years: 1993, 2013 and 2018 on 39 circular plots selected from a net of 400 m × 400 m plots regularly distributed within the whole area of the national park. The 13 plots were established in the subalpine mountain belt covered with Plagiothecio-Piceetum community and the remaining 26 were located in both lower and transition mountain belts within mixed spruce-fir-beech forests (Dentario glandulosae-Fagetum community) with dominating share of Picea abies in 1993. In 1993 ten spruces growing closest to the centre of each plot were selected and a survey of epiphytic lichens on a patch of bark was conducted from the ground level up to 2 m above the ground, around the trunk. The rate of coverage for each recorded lichen species on each patch was also estimated according to the following ranks: +: <1%; 1: 1%-5%; 2: 6%-25%; 3: 26%-50%; 4: 51%-75%; 5: 76%-100%. Because of the dynamic dieback of spruce forests within the Gorce range, mainly caused by the activity of European spruce bark beetle (Ips typographus L.), final analyses and interpretations were carried out using phytosociological data collected on 33 plots from 186 live spruces (out of 370 trees analysed in the first stage of research). Full species list with the sum of cover coefficients is available in Appendix 1. In the years 1993 and 2018 at each plot (wherever possible) thalli of Hypogymnia physodes were collected from the branches of standing trees or recently fallen branches or crowns of windthrowed trees (branches with needles still green, not in contact with the soil) and analysed for the concentrations of accumulated sulfur and heavy metals: Cr, Mn, Cu, Zn, Ni, Cd and Pb. In 2018, in addition to the previous set, concentration of nitrogen was also measured. The lack of suitable substrate in the immediate vicinity occurred in 3 of 33 sites in 2018. Eventually, the comparative analysis of the IAP values and the concentrations of accumulated elements were carried out for the remaining 30 sites.

Laboratory work
Lichen specimens, for which the identification of species based on morphological features was impossible in situ, were further investigated in a laboratory. The specimens were identified using standard microscopic techniques including simple spot test reactions with KOH, CaOCl and an alcohol solution of paraphenyldiamine (PPD). The nomenclature of lichen taxa followed Fałtynowicz and Kossowska (2018). Taxa were identified to the species level, except for sterile representatives of the genus Lepraria, which share would have been difficult to estimate because of morphological similarity. The dominant epiphytic representative of the genus in this study area was Lepraria jackii, but also L. elobata and L. finkii were frequent (Czarnota and Kukwa 2001). Currently recognized species in Cladonia pyxidata, Micarea micrococca and M. prasina complexes were also not precisely identified due to the unclear nomenclature or taxonomic status of cryptic species within these groups at the beginning of the study in 1993 (Wirth 1995;Smith et al. 2009). Finally, the taxa were treated as species sensu lato.
In 1993 the concentrations of heavy metals in Hypogymnia physodes thalli samples (ca. 1.2 g of dry weight) were analysed using Atomic Absorption Spectrometry (AAS) after preceding mineralization in 65% nitric acid. The measurement uncertainty was estimated in accordance to the upon standard procedure based on 7 replicates from one randomly chosen sample from site no. 21086. The relative standard deviation based on the series of measurements was: 3.7% for Cr, 1.59% for Mn, 2.83% for Cu, 0.64% for Zn, 32.62% for Ni, 5.03% for Cd and 5.06% for Pb. The concentration of sulfur was determined by nephelometry (see Czarnota 1995).
Concentrations of heavy metals in H. physodes thalli collected in 2018 were also analysed using methods described above. The averaged absorbance values of three technical repeats for each element were used to determine the concentrations, although a qualitative analysis of sulfur and nitrogen was performed using the CHNOS instrumental method. Samples of a homogeneous mixture of manually cleaned thalli collected from each site, weighing approximately 20 mg, were placed in tin foil containers and processed with the vario EL cube elemental analyser according to the procedure developed by Namieśnik and Jamrógiewicz (1998). Cleaned gaseous combustion products undergo separation in absorption columns into individual components (nitrogen, carbon dioxide, sulfur dioxide, water vapour) and are detected by the measuring cell of the TCD thermal conductivity detector. Each sample was divided into three subsamples which were measured thrice in order to obtain the value of standard uncertainty for each site. Both methods were calibrated and validated on the basis of certified reference material (SRM 1547 peach leaves). Recovery rates of N and S in the elemental analysis were 94% and 82%, respectively. In the AAS method, the recovery of Zn and Cu exceeded 90%, while for Mn, Pb, Cd, Ni and Cr >80% recovery was obtained.

Index of Atmospheric Purity
The phytosociological data from each study site were used for the calculation of the Index of Atmospheric Purity (IAP) in accordance to the original methodology of LeBlanc and De Sloover (1970) by using the following formula: where, n is the number of species; f is the frequency-coverage score of each species calculated by compiling the constancy class of occurrence and the rate of cover in the analysed plot of sampled trunks (see Table 1); Q is ecological sensitivity Table 1 The assessment of f coefficient used in IAP validation; criteria for the constancy occurrence class of species: I: 20%; II: 21%-40%; III: 41%-60%; IV: 61%-80%; V: 81%-100% (the criteria for determining the f coefficient have been specified here by authors).
Index of patch cover (%) Constancy class f coefficient value index of each species calculated based on the mean number of species coexisting on a site with given species.

Ecological indicator values
Each of the epiphytic lichen species was assigned an ecological indicator value proposed by Wirth (2010), describing habitat requirements on the 1-9 scale with relation to light (L), nutrients (N) and substrate reaction (R). The higher scale values correspond with increased either light intensity, concentration of nutrients or pH of the habitat. A small number of lichen species that have not been previously characterised by Wirth were described using a similar 5-degree scale proposed by Fabiszewski and Szczepańska (2010) developed for species common in Poland. Note that the value '1' on this scale refers to the first two levels on Wirth's scale and the adopted value in this case was '1.5'. Lichen species not mentioned in these two publications were assessed indirectly based on the indicator values for co-occurring species. The overall value for a given site was the mean value of indicators for all species found on a study site.

Total Index of Pollution
The concentrations of heavy metals in H. physodes thalli collected from the study sites were presented in a standardised format as a Total Index of Pollution (Grodzińska 1978) using the formula: = ∑ , for: = ( − ̅ )/ ̅ , where: observed concentration of i-th element in j-th sample; -mean concentration of i-th element for all samples.

Graphical interpolation and statistical analyses
Spatial distributions of IAP values, mean Wirth's ecological indicator values and accumulation of examined elements in the H. physodes thalli were presented on the map of Gorce National Park using the Inverse Distance Weighted (IDW) interpolation method in Quantum GIS v. 2.18.20 software. As sulfur concentration had been measured by using different methodologies throughout years of study, the sulfur pollution comparison was neglected and the spatial distributions of this element were only compared.
Data were processed using Statistica v. 13.1 and Canoco v. 5. software. Generalised linear model with Poisson distribution and log link function has been adjusted to predict: the impact on IAP indices, mean values of L, R and N coefficients of Wirth's ecological scale determined for the three years, i.e. 1993, 2013 and 2018, seven trace metals and sulphur measured in 1993 and 2018, and nitrogen measured in 2018. This model is appropriate for abnormally distributed response data based on Braun-Blanquet scale and effective under the assumption that the lichen's response increases approximately exponentially as the concentration of toxic substances grows. The maximum-likelihood approach in the model for the 1993 data was used (when the over dispersion phenomenon did not occur) and the quasilikelihood approach in two latter study years (when dispersion was >1.5). Selection of best predictors in the models for individual years was done by backward elimination based on Wald test. The nonparametric Mann-Whitney U test was applied to compare the ranked ecological indicator values for different forest communities.
The effect of environmental variables on the epiphytic lichen biota associated with spruce (estimated on the Wirth's ecological indicator values as well as the concentrations of S, N and heavy metals in H. physodes thalli) was also checked based on species abundances for each of the three research periods. Redundancy Analysis (RDA) was performed to assess the influence of that data in explaining the variance of the coefficient f (from the IAP equation) calculated for each lichen species from a study plot according to the criteria shown in Table 1. Analyses were performed on the basis of 999 unrestricted, random permutations and the significance of the model, and the ordination axis were tested using pseudo-F statistics based on the sum of all canonical eigenvalues (Miller 1975). The response data were previously centred by columns and square transformed in order to normalize their distribution and limit the effect of most dominant species. Most influencing environmental variables were sequentially added to the model using forward selection method. Only those species, for which the frequency on a site in a given period was higher than 5%, were analysed (Table 2). Species, for which the variances had not been sufficiently explained by the fitted model on the basis of their response curves, were removed in the last step.
Variation partitioning procedure (Borcard et al. 1992) was performed between two distinct groups of tested predictors, i.e. ecological indices and chromium, to identify the importance of its simple and conditional effects on the response variables.
Regarding heavy metals, a trend towards their increased deposition in H. physodes thalli has been observed. A standardised total index for heavy metals (Sj) was in the range of 2.777-5.536 in 1993 and 3.602-8.712 in 2018. Significant changes in the concentrations of 4 analysed elements were found between the examined years. The mean concentration of Cr reduced from 5.48 ppm to 0.88 ppm. On the other hand, there was an effective increase in the concentration of other elements: Mn (from 42.67 to 235.59 ppm), Ni (from 2.18 to 5.83 ppm) and Cd (from 1.07 to 1.59 ppm).

IAP
The values of IAP on study plots in 1993 were in the range of 14.5-56.2 ( ̅ = 32.4). After two decades, in 2013, IAP was in the range of 6.6-50.7 ( ̅ = 27.5). The fall in IAP during that period was caused by a decreasing rate of trunk coverage with epiphytic lichens from 46.0% to 28.3%. A significant (α=0.05) decrease in share of dominant species was observed for: Lecanora conizaeoides, Hypogymnia physodes, Xylopsora caradocensis and sensitive foliose species from the genera Usnea and Parmeliopsis as well as Platismatia glauca, Pseudevernia furfuracea and Nephromopsis chlorophylla. In 2018 the mean IAP index increased again reaching values between 9.2-62.2 ( ̅ = 36.9). The trend observed in 2013-2018 was associated both with an increase in the number (from 11 to 13 per plot on average) and the appearance frequency of crustose and shadetolerant species from the genera Chaenotheca and Micarea. The mean trunk coverage (30.4%) was comparable to the one in 2013. The list of IAP indicator values in the three research years is given in Appendix 3.

Relationships between IAP, accumulation of elements and ecological indicator values
The research showed no significant correlation (α = 0.05; Table 3) between the IAP values and concentrations of both sulfur and heavy metals in the same plots in 1993, although manganese (r s = 0.27) and nickel (r s = -0.27) had a moderate impact. The cumulative impact of examined metals (Sj index) on IAP values in the first term was weakly negative (r s = -0.10). In 2018 a stronger but still insignificant negative correlation between IAP and concentration of S (r s = -0.24) was observed, whereas no correlation of IAP with the content of N occured. In contrast to 1993, the correlation of IAP with Mn in 2018 reversed completely, although it was still not significant (r s = -0.28). In the same time period, a positive correlation of IAP with Cd increased considerably (to r s = 0.33) and Sj index changed to weakly positive (r s = 0.19).
In 2013 and 2018 IAP values was influenced with the type of forest community, being higher in subalpine spruce forests than in mixed forests at lower altitudes (Figures 1 & 2). The spatial distribution of deposited heavy metals (Sj index) for both years indicate the multidirectional transport of these pollutants, mainly determined by the course of river valleys, where metals are deposited usually on mountain slopes at lower altitudes. The air contaminated with S and N mainly flows from the north and north-east, i.e. from industrialized regions of Upper Silesia and Kraków, as indicated by the increased concentrations of these elements on the windward slopes of the highest hills (Figure 1).
Knowing the environmental requirements of lichens, Wirth (2010) developed an ecological scale in which individual species are assigned optimal values of physical parameters of the environment, just like in the Ellenberg scale created for vascular plants (Ellenberg 1974;Ellenberg et al. 1992). The light indicator (L) for lichens surveyed in Gorce NP was in the range of 3-8. Most species were assigned value 5 (46% observations) corresponding with semi-shaded habitats (>10% of relative Table 3 Spearman's ranked correlation coefficient (rs) values between the IAP value and the mean of ecological scale indicator (L -insolation index, R -substrate reaction index, N -nutrition index; according to Wirth (2010) and Fabiszewski and Szczepańska (2010) intensity of natural light). It was found, that L values differ between sites depending on the type of forest community. The optimum values of light assigned to species found in subalpine spruce forests (P-P) were higher than that in the Carpathian beech forest of lower belt (Dg-F) and had significantly different distributions (Table 4). The L indicator in each study year correlated with IAP at the level of r s > 0.3 (Table 3), but the correlation was statistically significant only in 2013. The comparison of spatial distributions of IAP values and L indicator values for 2013 revealed that higher availability of light on the plots corresponded with IAP values mainly in subalpine spruce forests, although this trend is not persistent, as is indicated by the drop in the correlation magnitude in 2018. The R indicator value (range of 1.5-5.5) was found to be 2 in most cases (44%) and pointed to the dominance of very acidic habitats with pH of 3.4-4.0 (Wirth 2010). In 1993 the R median was higher in the subalpine spruce forest (P-P), whereas the relation reversed in 2013 and 2018 (Table 4). The correlation between IAP and mean R distributions in 1993 was not significant (Table 3), but in two latter study years the correlation values between these indicators ranks expanded, i.e. r s2013 = -0.54 and r s2018 = -0.30 (Table 3), becoming  a relevant trend in the second term. The indicator value for nutrients (N) ranged from 1 (habitats poor in nutrients) to 5 (mesotrophic habitats), where the most frequent values were 2 and 3 corresponding with loweutrophicated or non-eutrophicated habitats (47% observations in total), respectively. Higher N values characterised Carpathian beech forest (Dg-F) in each study year (Table 4), but significant differences were found in 2013 and 2018. In subsequent years the N indicator ranks significantly correlated with IAP (Table 3), and considering all analysed factors, the eutrophication had the strongest negative effect on the lichen biota associated with spruce. Increased IAP values corresponded to decreased mean values of N indices on a study plot and simultaneously the occurrence of more nitrophilous species correlated significantly negatively with cadmium contamination (r s = -0.38).

Generalized linear model for assessing accumulation of elements and ecological indicator values as predictors of IAP
For the further consideration of the effects of the individual predictors tested on the IAP values in each study year, generalized linear models were built based on the Poisson distribution (Table 5). Models fit were tested by using Pearson χ 2 test and analyses of residuals. The model for 1993 was the best fitted one (χ 2 /df = 1.07), where both L and N indices as well as cumulated nickel content were significant predictors of IAP (α=0.05). In the model for 2013 N index was significant wherein χ 2 /df = 1.82. Finally, in 2018 N index, Cd and Mn contents were relevant factors, whereas χ 2 /df = 1.64. Each time the N index was the most negatively related among other factors, indicating that nitrophilous species decreased their share in lichen communities during the last 25 years. At the same time, the reaction coefficient R was not influential in any of the models, and the light index L was significant only in the first research term. That is incompatible with Spearman's correlation of these variables, shown in Table 3, where IAP correlated positively and significantly with an increase of both photophilous (2013) and acidophilous (2013,2018) species. Accumulated nitrogen in 1993 turned out to be a better predictor for IAP values than photophilous species, despite being less correlated to it (see Table 3). The same was true for Mn and Cd, which abundances were more powerful in the model than slightly better correlated L and R indicators. The differences are resulting from an interaction effect occurring between predictors included in the model.

Multivariate analysis
RDA for individual study years identified significant environmental factors (α = 0.05), their directions and magnitude of effect on the presence of epiphytic lichens on spruce trees ( Figure 3A-C). The analysis of data for 1993 (L, R, N, sulfur, 7 metals, Sj) explained 43.0% of the total variance for species composition (34.8% of adjusted variance, p = 0.001; Tables 6 and 7). Factors with a significant influence included the three indicator values and chromium (Table 7). L and N indices were the most important factors, the first of which was negatively and the second positively correlated with the first ordination axis. R index was better positively related to second axis, and Cr accumulation was about equally related to both of them. Based on variation partitioning procedure, only the conditional effects of Cr were significant, but not their simple effects, meaning that they were not significant as an independent variable. The weak influence of chromium is also shown on the species-environmental variables biplot on the Figure 3A. The 2013 analysis tested ecological indicator values only (the concentrations of elements were not measured that year). The model was significant at the level of p = 0.001. All three indicator values had a significant effect on the species composition at that time, explaining in total 34.1% of data variance ( Table 6). The best correlation with the first two axis (Table 7) was found for reaction index (R), which is thus constituted the best predictor. L was highly correlated with all three significant axis whereas N with the first and the third of them. Among thirteen analysed variables in 2018 (L, N, R, sulfur, nitrogen, 7 metals, Sj), two were significant and explained together 37.8% of the total (33.2% adjusted) variance of the community structure (p = 0.001; Tables 6 & 7). These were L and N values, well correlated with two first RDA axis (Table 7). The concentrations of heavy metals measured either separately or jointly in lichen thalli in two study years, with the exception of Cr in 1993, had no significant effect on the frequency and constancy of lichen species.  Light ( For the purpose of air quality assessment in the studied sites, the IAP index was determined, based on the phytosociological parameters of epiphyte communities, such as: species composition, coverage of the examined bark patch by each species and the constancy of their occurrence. Despite the obvious weakness of treating all species equally, it is widely recognized as being effective in its role when it is known that they differ in tolerance to pollution and microclimatic conditions. The assumptions for its use are met, as long as pollution is a limiting factor for the occurrence of lichens. Otherwise, as in many ecosystems of natural forests away from emission sources, the determinant are microclimatic factors, modified by the continuity and diversity of the ecosystem and the dynamics of forest communities, as well as the substrate properties.

Discussion
The mean concentrations of S and N in Gorce range measured in lichen thalli in 2018 were comparable to those measured with the same methods in 2015 in samples of the same lichen species transplanted to the area of Rzeszów (a city of a population of 200.000), located ca. 200 km east of Gorce range, reaching 0.133% d.w. and 1.248% d.w. of sulfur and nitrogen, respectively (Tanona et al. 2017). Similar results were reported that time in a control sample collected from large Solska Forest, i.e. 0.144% d.w. and 1.22% d.w. of S and N in thalli, respectively.
In 2018 the concentrations of trace metals (except Ni) in lichen thalli from the GNP were slightly lower than the levels measured by Klimek et al. (2015) in thalli of H. physodes collected in five towns located in the Western Beskidy Mts. The heavy metals pollution level in Gorce National Park in both periods was lower than in nearby Kraków (Białońska and Dayan 2005) and Kielce, located ca. 200 km N of Gorce range (Jóźwiak 2007), but higher than in Rzeszów (Tanona et al. 2017). The above results suggest that the inflow of pollutants from Upper Silesia and Kraków regions might have an important effect on the air quality in this range.
Considering the well examined, clear and adequate biological response of lichens to the presence of toxic elements and compounds (SO 2 , NO x and heavy metals) in urban areas (e.g. Ranta 2001;Bačkor et al. 2003), the zones designated by the interpolation of IAP values in the forested area of Gorce National Park could, at least partially, correspond with the accumulation of total sulfur, nitrogen and/or the content of metals in the sampled thalli. In the light of our results, i.e. the lack of clear correlation between IAP and air pollutants accumulated in lichen thalli, the use of this method for the assessment of air pollution in large mountain natural forests seems to be at least controversial.
At the same time, the phytosociological IAP index was better correlated to Wirth's ecological scale indicators. This most likely results from the positive response to forest decomposition caused by natural disasters in recent years, such as bark beetle outbreaks and hurricane winds.
The important effect of light on the biodiversity and distribution of lichens in forests, which has been recently emphasized i.a. by Sales et al. (2016), Sevgi et al. (2019) and Vondrák et al. (2019), showed up here in the similarity between spatial distributions of IAP values and L indicator values for 2013.
An increase in correlation between IAP and R index was observed from the first towards the second and third study year, while the mean concentration of sulfur measured in H. physodes thalli was decreasing. In 1993 S content was almost two-fold higher than the concentration measured in 2018. This shows the simultaneous influence of other factors modifying the response of lichens growing on the sites and explains the lack of correlation between IAP and mean R distributions in 1993 in the first study year.
As it was previously said, the eutrophication (N index) had the most harmful effect on the studied lichen biota (Table 3). The negative correlation between both N and Sj indices and Cd concentration, with the simultaneous lack of a relationship between the N index and nitrogen deposition, indicate that the increase of IAP values was caused by the presence of species tolerating higher concentrations of heavy metals, particularly Cd, and not by the share of nitrophilous species, although this relationship is not statistically significant (r s = 0.19; Table 3). The lack of a significant relationship between the IAP indices and the concentrations of the examined metals in lichen thalli on the same forest plots show that the long-term (25 years) dynamics of the frequency and abundance of lichen species did not significantly depend on heavy metal contamination. It means that the level of these elements was and still lies within the tolerance of species inhabiting the monitored forests. Generalized linear model reveals that N index was the best predictor for IAP in three study years (Table 5). Although it highlighted such metals as Ni, Mn and Cd, the impact of accumulated elements on the IAP index is still uncertain and difficult to interpret due to the instability of the direction of interaction (e.g. Mn, see Table 3) as well as surprisingly positive relationships (Cd; see Table 3).
The content of S and N in sampled H. physodes thalli did not correspond with ecological indicator values describing the share of acidophilous and nitrophilous species (Table 3). The amount of pollutants absorbed and retained by lichens depends on various factors, including: the lichen species, the type of pollutant as well as climate and environmental conditions (Sett and Kundu 2016). Moreover, the deposition of airborne particles on the tree bark surface may depend for example on: microclimatic conditions, topography (altitude, exposure, slope gradient) or tree stand density (Kowalkowski 2003;Xu et al. 2016). The pH of bark, apart from its natural value, is not determined exclusively by the level of sulfur compounds adsorbed on the bark surface. It can also be reduced by nitrogen compounds, including ammonia and ammonium ions deriving from decomposed organic matter and carbonic acid being the product of the reaction between H 2 O and CO 2 (Brimblecombe 1994). Consequently, the spatial distribution of the elements accumulated in lichen thalli may differ from the content of these elements in the air (Liu et al. 2016;Węgrzyn et al. 2016). It shows that the concentrations of toxic elements in sampled lichen thalli could not be used as an ecological indicator in all environmental circumstances, especially in large forests.
Multivariate analyses demonstrated that sulfur and nitrogen deposited in lichen thalli, at levels measured in 1993 and 2018, have no relevant effect on the phytosociological parameters of lichen epiphytes associated with spruce and among heavy metals, only chromium had a certain impact. In each research term the availability of light and nutrients in a given habitat had a substantial effect. Lichen community structure was determined by availability of light in the first and the third study year and by substrate reaction in 2013. Comparison of these results with the generalized regression model related to the IAP indicates that light and reaction are responsible in particular for the species composition of the plot, while nutrients influence rather the diversity and abundance of species.
The IAP value decreased significantly over first 20 years (1993-2013) of observation. During the last period 2013-2018 the IAP value increased again exceeding its initial mean value (see Figures 1  & 2). A final explanation for the frequency changes and the richness of lichen species associated with Norway spruce should not be sought in the dynamic of air pollution (see Table 3) but in the dynamic of tree stands that took place over the last 25 years in the area of Gorce NP. These changes are mainly driven by bark beetle outbreaks and windstorms. Light availability in the forest may be decisive for the toxic effect of nitrogen at low concentrations (Carter et al. 2017) as well as influence the negative correlation between IAP indices and mean values of L and N indicators. Dynamic changes in the lower belt forests of the Gorce range reduce the availability of light on the forest floor because of the vigorous development of undergrowth. It is the main cause for the natural evolution of species composition towards communities tolerant to: shading, higher doses of nutrients and lower acidity of the substrate, which is reflected in higher ecological indicator values N and R. These conditions are suitable for lichen species, which throughout the study period were characterised by the highest increase of frequency in Dg-F forests: Coenogonium pineti (L = 3, R = 4, N = 4), Micarea micrococca (L = 4, R = 1.5, N = 3.5) and Chaenotheca stemonea (L = 3, R = 3, N = 2). In subalpine spruce forests P-P gaps in post-bark beetle stands persist much longer and pollutants can be leached by rain, which is more frequent in this zone. These conditions promote the expansion of oligotrophic, light-demanding species that prefer acidic phorophyte, for example Hypocenomyce scalaris (L = 6, R = 2, N = 2), Lecidea nylanderi (L = 5, R = 2, N = 1) and Violella fucata (L = 5, R = 3, N = 3). A significant decrease in the share of foliose and fruticose macrolichens in subsequent study years was observed mainly in the Dg-F zone and it can be largely attributed to the gradual shading of study sites. The decrease in the share of these species (less significant) in the P-P forest community was most likely caused by sudden changes in microclimate (insolation and humidity) following natural disasters that affected the surviving trees. A more frequent exposure of thalli to direct insolation causes greater evaporation of water from the bark, and thus, as reported by Gauslaa (2014), may consequently limit the growth of species that require higher humidity.
Many researchers have emphasized that the dynamics of lichen biota are driven by many environmental factors (including the most important ones, such as climate, microclimate and substrate properties (e.g. Asta et al. 2002;Giordani 2007;Agnan et al. 2017) as well as human pressure (including transformation of habitats; Chuquimarca et al. 2019), which should not be ignored while interpreting data on the diversity of lichens. At the same time, it should not be forgotten that factors of different origin interact in the environment, so the influence of "natural" factors may also be modified by the anthropogenic factors, even of far-reaching ones. Asta et al. (2002) reviewed methods for calculating phytosociological lichen indicators and proposed a modified methodology for the assessment of the lichen diversity value, namely LDV. According to the authors, this methodology has a wide range of applications, including nonurban areas, and allows for the long-term impact assessment of environmental stressors, including air pollution, eutrophication, anthropogenic environmental changes and climate change. Agnan et al. (2017) researched natural forests in France and Switzerland and reported only a weak negative correlation between the level of accumulated metallic elements and the LDV index, but no such relationship referred to IAP values, although IAP, LDV and abundance of lichens correlated better with the ecological indicator values for reaction, eutrophication, continentalism and light. Käffer and de Azevedo Martins (2014) used the original methodology of IAP while investigating virgin-type forests of Brazil and concluded that it was an adequate tool for the identification of areas valuable in terms of lichen biota.
The hereby results correspond well to the above studies, confirming that the original version of the IAP bioindication method can be applied to assess ecological conditions in large Carpathian forests, however not in their validation focused on environmental pollution. A detailed analysis of a range of environmental variables, that could potentially influence the examined epiphytic lichen biota, was not the aim of the present study, but will be the subject of another research paper by authors.
The strengths of using the IAP index in assessing the quality of ecological conditions in forests with a low level of air pollution are: i) the possibility of designating zones of high ecological value in forests, i.e. species rich, weakly affected by human, with preserved continuity of ecological processes; ii) the possibility of assessing the impact of stand dynamic, including various types of disturbances -natural and anthropogenic, on a group of lichens known to be sensitive to changes in their habitats; At the same time, there is a weakness of the method, namely the fact, that the value of the indicator may be determined by tolerant species, occurring in high abundance and with a high constancy. Nevertheless, its equation was constructed in such a way that it was more sensitive to the species diversity of the site than to changes in the surface coverage.

Conclusions
IAP method is recommended to use for the assessment of ecological conditions in forests. This particularly refers to large natural Carpathian forest, where the concentrations of toxic compounds and elements (to a certain level) are assumed not being the most important criterion determining the growth of lichens. IAP is suitable for the identification of forest patches most valuable in terms of the protection of epiphytes but also for the long-term monitoring of the response of epiphytic lichen biota to local environmental changes induced both by humans and natural causes. with significant logistical assistance of the Gorce National Park. Open access funding provided by University of Rzeszów, within the CRUI-CARE Agreement. Special thanks go to Mateusz Czarnota for proof reading.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.