The multi-isotope biogeochemistry (S, C, N and Pb) of Hypogymnia physodes lichens: air quality approach in the Świętokrzyski National Park, Poland

ABSTRACT The isotope biogeochemistry of bioindicators has widely demonstrated its added value in environmental issues by allowing to precisely identify sources of contamination. Most of the studies are based on studying one or two isotope systematics. Here, we are presenting an innovative multi-proxy approach that combines chemistry with both stable (C, S, N) and radiogenic (Pb) isotope systematics. Using Hypogymnia physodes bioindicators, we evaluated air quality in the complex environment of the Świętokrzyski National Park (ŚNP, Poland) with the ultimate objective of isotopically identifying the sources responsible for the observed contamination. Combining the isotope systematics showed that home heating is a major source of contamination in winter, whereas the contribution of road traffic increases during the summer. Pb isotope ratios identified industrial activities as the major source of this metal in the atmosphere.


Introduction
Anthropogenic activities are generating contaminants and pollutants that have now become a major issue when trying to preserve our environmental quality, and these can be studied using bioindicators [1,2]. Environmental contaminants affect organisms such as plants, which at first lead to bio-accumulation of these contaminants and eventually to changes in their anatomy and physiology [3]. Thus, the physiological and morphological impacts of air pollutants on living organisms have been widely investigated for the past decades [2,[4][5][6][7][8][9][10]. The organisms that possess the ability to record pollution events are referred as bioindicators. They present many advantages and in particular their potential to bio-concentrate both organic and inorganic contaminants, which makes possible the use of bioindicators to assess air quality [6].
Lichens represent a large group of reliable bioindicators, due to their ability to accumulate chemical elements in amounts usually exceeding their physiological needs [e.g . 11]. They are thus able to survive by maximising precipitation and atmospheric aerosol uptake [12]. Lichens are organisms of a dual nature that are made up of two different componentsphotobionts and mycobionts. The interaction of these components, their growth, nutrition and reproduction is extremely complex and depends on many environmental factors [13], which makes it difficult to thoroughly study and understand the physiology of lichens, and thus the processes of accumulation of various compounds in their thallus. It is known that nitrogen and sulfur compounds, heavy metals etc. present in lichens are absorbed by the thalli to different degrees, also depending on the level of environmental pollution. These compounds are not permanently accumulated in the lichen thallus, and their concentrations do not increase proportionally with time. The absorption of some compounds may be inhibited by specific lichen substances that enable its physiological processes [14]. The content of substances absorbed by lichens in their thallus may change within the annual cycle, depending on various factors such as temperature, sunlight, oxygenation, humidity, as well as elevation [15,16]. Lichens are physiologically most active during episodes of high humidity (rainfall, fog, dew) and when their thallus grows. The activity of lichens ceases in dry periods as poor hydration of the thallus inhibits the process of photosynthesis in algae [17][18][19]. After intense physiological activity during wet months, the lichen turns into a state of limited activity due to too low temperatures (below 0°C) or too high temperatures (causing the drying of the lichen thallus). Still, one of the major limitations of this bioindicator approach is that it yields information that is time-averaged over the lifetime of the bioindicator, which in turn is often difficult to determine [20]. Lichens are often used for determining the levels of pollutants and for tracking their corresponding sources: for example, sulfur (e.g. [21,22]), nitrogen (e.g. [23][24][25]), heavy metals (e.g. [6,26,27]).
With that in mind, we are here proposing to use an innovative multi-proxy approach, coupling chemistry with both stable (C, N, S) and radiogenic (Pb) isotopes from bioindicators (Hypogymnia physodes), to characterise air quality in a complex environment. To fulfil this objective, we selected the Świętokrzyski National Park, Poland, a supposedly rural study site that has been proved to be also impacted by short and long-term anthropogenic contamination [10]. We specifically aimed at answering the following questions: (i) What are the source(s) and process(es) controlling local air quality? (ii) Is there any seasonality?, and (iii) Can this multi-proxy approach be up-scaled and implemented elsewhere?

Study area
The Świętokrzyski National Park (ŚNP; Figure 1) is located in the southeastern part of Poland, within the Świętokrzyskie Voivodship. Many sources of air pollution are present in the Voivodship. The most significant ones are the industrial sector, road traffic and home heating for local single-and multi-family buildings. The air pollution in the study area might relate to the metallurgical industry as steel plant, and zinc smelter are located nearby. Also, sulfur chemical plants can contribute to this situation, as well as coal-fired heat and power stations. Air quality in the study area is also impacted by air pollution originating from the neighbouring Voivodships (mainly from the Silesian and the Lesser Poland Voivodships [42]). The main nearby sources of air contamination are located within the west and south parts of the Voivodship (Figure 1). Within the study area the prevailing winds blow from the south (45%) and the west (40.2%) while the southeast and east directions contribute to 8.3 and 4.8%, respectively [43].
The meteorological station, located in ŚNP, reports a mean annual temperature of 5.7°C , compared to 7.7°C for the whole Świętokrzyskie Voivodship. The warmest month is July and the coldest one is January, with mean temperatures of 15.7 and −4.6°C, respectively [44]. Between 1961 and 1990, minimal precipitations we generally observed in winter (February) and the maximal ones during the summer (generally in July; [44]). The annual sum of precipitations reaches 800 mm in the Świętokrzyskie mountain range, and 550 mm in the Nidziańska Basin. The highest mean sum of precipitations (823 mm) is observed at Święty Krzyż [44].

Materials and methods
3.1. Sampling strategy 20 sampling sites were selected (Figure 1), supposedly reflecting distinct environmental conditions. Sites 9, 10 and 19 (located close to roads) were expected to be mostly impacted by local road traffic, 11, 12 and 15 by long-range pollutants, and 3, 4, 6, 7 by home heating. Site 8 was supposed to represent a background site as it is located in the centre of the national park with no roads nearby. The rest of sites was considered to be influenced by various emission sources, including industrial activities, road traffic or home heating. The altitudes of the sampling sites varied between 284 and 597 m a.s.l. Three of them were higher than 450 m a.s.l.: sites 11 (484 m a.s.l.), 12 (597 m a.s.l.) and 15 (591 m a.s.l.) ( Figure 1).

Bioindicators
We elected to use Hypogymnia physodes (L.), an epiphytic lichen as this species has widely been used as a bioindicator for air pollution (e.g. [10,26,[45][46][47]). Lichen samples were collected during two campaigns: winter (1-3 February 2013) and summer (11-13 June 2013) inside or in the close vicinity of the National Park ( Figure 1). Between 1 and 20 g of Hypogymnia physodes were collected at each sampling site from the trunks and branches of Abies alba trees, at a height of 1-2 m, using a plastic knife. We selected similar sizes of branches which were supposed to provide the lichens of similar age. For each sampling point we have sampled a couple of different trees, and the overall collected material was homogenised to obtain representative lichen samples. After collection, the samples were packed in paper bags and then transported to our laboratory where they were weighted.

Analytical methods
Lichen samples were homogenised, and an aliquot of ∼6 mg for nitrogen and of ∼2 mg for carbon were inserted into two separate tin capsules that were then sealed for their corresponding concentration and stable isotope analysis. Carbon and nitrogen concentrations were determined using a Carlo Erba NC 2500 elemental analyzer. The corresponding carbon (δ 13 C) and nitrogen (δ 15 N) stable isotope compositions were measured using a Vario Micro CubeTM elemental analyzer coupled to an Isoprime 100TM isotope ratio mass spectrometer (IRMS). The error on the δ 13 C analysis was <0.2 ‰ and of <0.2 ‰ for δ 15 N. For sulfur (S adsorbed ), lichen samples were washed with deionised water for one hour. The obtained solution was filtrated, heated at 100°C and acidified down to pH = 2 with an 18% HCl acid solution. An excess of BaCl 2 was then immediately added to precipitate S as BaSO 4 . The formed BaSO 4 was purified by multiple rinsing with deionised water in order to remove remaining chlorides. Finally, BaSO 4 was collected after drying at 50°C. Initial sulfur concentration in the samples were calculated from the generated BaSO 4 . For measuring the corresponding δ 34 S, BaSO 4 was totally converted into SO 2 during reaction with V 2 O 5 and after that, it was cryogenically purified off-line. We followed the procedure described in [48]. Finally, δ 34 S was determined on the pure SO 2 gas using a Delta Plus Advantage IRMS. The results of δ 34 S were expressed against the international Vienna Cañon Diablo Troilite (VCDT) scale. Isotope data were normalised using international IAEA standards i.e. NBS-127, IAEA-SO-5and IAEA-SO-6. The error on the δ 34 S analysis was <0.2 ‰. Organic sulfur in the lichen samples was pre-concentrated using the Parr Bomb (1108 Oxygen combustion vessel of a Parr Instrument Company, USA) technique, where the sample is mineralised by combustion under high pressure (2.5 MPa) of oxygen. The solution obtained after the mineralisation was then filtered, heated and acidified using 18% HCl down to pH = 2. BaCl 2 was then added in excess to convert dissolved sulfates into solid BaSO 4 . In order to remove formed chlorides, BaSO 4 were purified by multiple rinsing with deionised water. Samples were filtered, and the collected BaSO 4 was burnt in porcelain vessel at 800°C for approx. 4 h (following the internal procedure of the Thünen Institute in Braunschweig). %S org was calculated from the obtained BaSO 4 , from which δ 34 S was determined. For lead lichen samples were first aciddigested and chemically purified following [49], through ion exchange (AG1-X8 resin) using ultra-pure HBr/HCl. The Pb isotope ratios ( 206 Pb/ 204 Pb, 207 Pb/ 204 Pb, 208 Pb/ 204 Pb, 207 Pb/ 206 Pb and 208 Pb/ 206 Pb) were measured by multicollection inductively coupled plasma mass spectrometry (MC-ICP-MS) on a NuPlasma II with an Aridus II desolving membrane as the introduction system. All samples were doped with thallium to correct for the mass bias (e.g. [50]). Measurement of NBS-981 every five samples throughout the analytical sessions allows to get a Tl-Pb mass bias relationship [51]. The measurements of the in-house CGPb-1 standard yielded a long-term average (+/-2σ, n = 22) 208 Pb/ 206 Pb = 2.0425 ± 0.0003; 206 Pb/ 204 Pb = 18.596 ± 0.008, 207 Pb/ 204 Pb = 15.702 ± 0.005, 208 Pb/ 204 Pb = 37.98 ± 0.02, 207 Pb/ 206 Pb = 0.8444 ± 0.0002.

Statistical analysis
A first conventional correlation analysis did not yield satisfying, easy interpretative correlations, so we opted to perform a principal component analysis (PCA) on our sample dataset. All calculations were carried out using the Statistica 13 software [52]. As we expected lichen major elements to be controlled by both natural and anthropogenic factors/processes, the PCA approach was selected to help discriminate them. Manly [53] demonstrated that PCA requires a limited number of useful variables and helps describe relations amongst them. All parameters (i.e. height, CO 2 , NO 2 and SO 2 levels, C, N and S contents, as well as their respective stable isotope compositions) displayed normal or close-to-normal distributions (Shapiro-Wilk test applied, [54]). This demonstrated the representativeness of our data and yielded accurate results in the PCA [53,55]. For winter and summer samples, five factors and their shares (scree test), Factor loadings (varimax rotation) and Factor scores were calculated [55][56][57].

Results
Results for both winter and summer campaigns are available in the Supplementary Material section.

Lead isotope ratios
During winter 206

Organic carbon, nitrogen and sulfur, and adsorbed sulfur in lichen samples
Nitrogen concentrations in the lichens, %N org , ranged from 1.2 to 2.6% (average of 1.9 ± 0.5%) for summer, and from 0.9 to 3.1% (average of 1.9 ± 0.7%) for winter. %C org varied from 35.1 to 47.0% (average of 43.5 ± 2.6%) in winter, and from 43.3 to 49.4% (average of 46.4 ± 1.7%) in summer. As lichens are recognised as sensitive sulfur biomonitors [6,58], we determined the levels of both inbuilt sulfur (S org ) and S deposited on the lichen surface (%S adsorbed ). S org varied between 690 and 2060 ppm (average of 1180 ± 350 ppm) in summer (Table S5), and between 60 and 1010 ppm, (average of 450 ± 270 ppm) in winter (Table S4). On the other hand, S adsorbed ranged from 130 to 1530 ppm (average of 476.5 ± 320 ppm) in summer (Table S5) and from 60 to 1370 ppm (average of 406.5 ± 320 ppm) in winter (Table S4).

Nitrogen
Previous studies reported an overall range of nitrogen concentrations in lichens from 0.32 to 4.69% for different lichens species collected worldwide [2,8,[23][24][25]45,[59][60][61][62]. For our samples, %N org ranged from 0.9 to 3.1% in winter (average of 1.9 ± 0.7%) and from 1.2 to 2.6% in summer (average of 1.9 ± 0.5%). The absence of visible seasonality in the %N org probably resulted from the very low yearly growth rates of lichens, which are known to yield N concentrations averaged over several years [23]. However, individual samples show that some lichens display similar %N org values for both seasons (6,8,15,16,18,20), while for others slight seasonal differences can be detected (Figure 2(A)). This may be explained by: (a) seasonal variations in the number and types of contamination sources, or (b) seasonal variations in the lichen assimilation rates. Interestingly, the highest amounts of assimilated N org , for both seasons, were observed for lichens collected at the highest altitudes (11,12, and 15green square; Figure 2(A)). This can be explained by specific meteorological conditions (higher precipitations and moisture content) that enhance organic nitrogen (N org ) assimilation and accumulation, even during winter. δ 15 N isotope compositions of the major N atmospheric contaminants are usually representative of their corresponding emission source. It has been widely demonstrated that δ 15 N org of lichens are strongly linked to that of the corresponding emission source(s) and are N speciation-dependent. Here, the δ 15 N org of lichens were consistent with the range previously reported in diversified lichen species (e.g. −16.0 to +4.3 ‰; Table S1; [2,24,25,30,45,63,64]). The main N-compounds controlling the δ 15 N of lichens are: NO x from road traffic, ranging between −28.1 and +17 ‰ (e.g. [65][66][67][68][69][70]), ammonia from road traffic (δ 15 N ranging from −5 to −3 ‰ [41]) or N-bearing dust deposited on surrounding roads (δ 15 N≅+6.8 ‰ [67]). Hence, lower atmospheric δ 15 N values are usually observed in areas dominated by agriculture, especially in presence of fertiliser and manure spreading [71][72][73] whereas higher values are associated with industrial activities and road traffic in urban areas (e.g. [67,71,74]). In our study, samples collected close to roads (9, 10 and 19) showed a slight increase in their %N org and δ 15 N org in summer compared to winter (Figure 2; purple squares). This came from higher NO x emissions from summer road traffic due to a higher tourist attendance during this season in the Świętokrzyskie area [75]. Other samples (e.g. 3, 4, 6, 7), presenting lower δ 15 N org for summer when comparing to winter (Figure 2(B); orange square) might be connected with an increased application of fertilisers in local agricultural fields (especially in April; [76]) and increased precipitations, which enhance lichens' physiological processes during that period. Additionally, this may be enhanced by the high concentrations of ammonium in rainfall in Poland in April/May [77]. On the other hand, as lichens are located close to inhabited area, winter samples may reflect the early start of the heating season, leading to higher δ 15 N org . The highest differences between seasons were found at sampling locations 19 (intense road traffic in summer) and 6 (home-heating in winter and agriculture impact

ISOTOPES IN ENVIRONMENTAL AND HEALTH STUDIES
before summer). The rest of the sampling locations, including those located on the equity line in Figure 2(B), can be attributed to the mixing of several sources, which are difficult to identify at this stage.

Carbon
Carbon contents in the lichen samples differed between the winter and summer seasons, with lower %C org values in winter (Tables S4 and S5). Spatial variations were also observed with higher concentrations along roads and on peaks during the summer season, consistent with the findings of e.g. [78] and [8]. The values of %C org in lichens can generally range from 8.9 to 49.4% [2,8,45,[78][79][80][81], reflecting the influence of several factors, both natural (light, moistness/humidity, CO 2 derived from root respiration) and anthropogenic (CO 2 from fossil fuel combustion).
The average δ 13 C org of lichens was −26.6 ± 0.7 ‰ for winter and −26.9 ± 0.9 ‰ for summer. Isotopic compositions in lichens are influenced by the type of sources emitting atmospheric carbon and their corresponding δ 13 C ranges (Table S3). Apart from biogenic and anthropogenic emissions the final δ 13 C will be also shaped by other parameters such as local ozone concentrations, elevation, humidity and rainfall [29]. As shown in previous studies, increased rainfall leads to a depletion in 13 C and consequently lower δ 13 C org [28,29,82]. Ozone is present in higher amounts at higher elevations, and it is known to cause the decrease of photosynthesis [29], leading to the 13 C enrichment (e.g. [29,45,83]). However, we did not observe such a correlation between δ 13 C org and the elevation in our study (R 2 = 0.07 in winter and R 2 = 0.01 in summer). The lack of correlation may be linked to the relatively low differences in elevations among our sampling sites (maximum of ∼300 m). Besides, the canopy effect may also have buffered our δ 13 C values as it has been shown to lead to 13 Cdepleted isotopic compositions for samples collected near ground surface [84]. Unexpectedly, the values of δ 13 C org for lichen samples collected near roads (samples 9, 10, and 19) did not clearly reflect the typical carbon isotope signal of fossil fuel combustion. This may be explained by the fact that they were collected within the ŚNP, where natural factors and processes may isotopically overprint the influence of nonpoint sources such as road traffic for carbon.
In this area in 2013, SO 2 concentrations were a few times higher in the winter than in summer [10], and as SO 2 represents a major source of sulfur for lichens (e.g. [5,21,87,88]) higher %S org values were hence expected in winter samples. However, almost all sampling points in winter were characterised by much lower %S org than in summer (Figure 3(A)), indicating that there must be a factor inhibiting the SO 2 assimilation by lichens in winter. We hypothesise that lichen sensitivity to SO 2 plays a major role in the assimilation of S. The hypothesis is supported by the fact that the dissolution of

ISOTOPES IN ENVIRONMENTAL AND HEALTH STUDIES
atmospheric SO 2 in rainwater produces acidic ions that are readily absorbed through the lichen thalli, disrupting the photosynthesis process [89].
Another parameter seemingly discriminating our results is the sampling location. We observed lower %S values for samples taken in the centre of the ŚNP (e.g. points 8 in summer, 1-7 and 10 in winter) and higher ones in its southern area (e.g. points 11, 13, 15 in summer, 12 in winter). This can be explained by the fact that (i) most of the emitting sources are located on the south and western parts of our study area (Figure 1), and (ii) the prevailing winds blow from S and W. This indicates that the ŚNP, as a mountain range, stops air masses, which results in higher %S values in its southern parts, especially at stations located on peaks (e.g. 11; Figure 3(A)), where the long-distance pollutants have the biggest influence. A high value was also found for point number 20, located outside the park with no natural barrier, which might also indicate the impact of long-distance pollutants.
The assimilation of S compounds can be assessed by the comparison between the δ 34 S org (inbuilt in lichen) and the δ 34 S adsorbed (adsorbed on lichen thalli). Here, seasonal differences exist for those two parameters (Figure 3(B)). Generally, δ 34 S of anthropogenic emissions are centred around +5 ‰, but still depend on the type and nature of the emission source (Table S2). In our study area, the possible S inputs include emissions from anthropogenic sources and precipitations, while other natural sources such as volcanic gases, marine aerosols, sea spray may be considered negligible. Hence, the δ 34 S values of these latest are not presented in Table S2. Winter samples were enriched in 34 S compared to those from summer, and the same trend is observed for the SO 4 2adsorbed on lichens. This may result of intensified coal combustion and/or industry processes during winter, leading to higher δ 34 S org values. It should be noted that isotope fractionation happens during S-compounds assimilation. In almost all samples the S isotope composition of the lichen itself is depleted in 34 S compared to the corresponding isotopic composition of the SO 4 2adsorbed on it. This might indicate a preferential assimilation of lighter isotopes by lichens that requires less energy, similar to what has been previously reported for other plants [90,91]. Interestingly, while we observed that %S org contents in lichens thalli in winter were much lower than in summer, there is a visible shift of δ 34 S org in winter towards more positive values. It suggests that despite high atmospheric SO 2 concentrations in winter lichens assimilate less S org due to their sensitivity to S compounds, but we still observe a shift in their δ 34 S related to the winter-specific sources of pollution. Moreover, in winter δ 34 S values yield a narrower range compared to summer ( Figure 4) and that can also be linked to the intensification of coal combustion in winter, resulting in higher δ 34 S. The larger range of isotopic compositions in summer hints at a possible ternary mixing between (i) dissolved SO 4 2in precipitations, and SO 2 from exhaust gases (ii) before and (iii) after desulfurisation ( Figure 4). Additionally, the maximum of precipitations was observed in July. Hence, the lichen isotope compositions in summer should be consistent with those of dissolved sulfate in rain.

Lead
Pb concentrations in lichens roughly increase with altitude ( Figure 5(A), R 2 = 0.57). These specific locations, as mentioned above, are especially subject to the long-range transport of pollution from surrounding industrial areas. Additionally, the amount of precipitations increases with the altitude, inherently increasing moisture in lichens, which ultimately enhances their heavy metal uptake [92]. However, we did not observe a positive relationship between altitude and Pb isotope ratios in lichens, contrary to the findings of Doucet and Carignan [37]. Similarly to the sulfur results, higher Pb concentrations were observed in the southern part of the park, again related to the ŚNP mountain range in its central part acting as a natural barrier. This is in agreement with the dominating winds in this area (S and W) that cause an inflow of air masses from the nearby industrial areas ( Figure 1). Hence, industrial activities are considered the major source of atmospheric Pb in the study area. This hypothesis is also confirmed by the average Pb isotope ratios in our lichens that is consistent with previous values reported for a French industrial zone. This comforts our conclusion that the ŚNP is impacted by industrial Pb. Additionally, Sensuła et al. [93] reported a 206 Pb/ 207 Pb range between 1.15 and 1.17 for pine needles in a highly industrialised region of southern Poland, consistent with our results.
The comparison of our lichen 206 Pb/ 207 Pb and 208 Pb/ 206 Pb ratios with the potential sources of atmospheric Pb shows that our results are encompassed within the typical isotope ranges for Pb in soils, minerals (e.g. galena, limestone), gasoline, lignite, coal and particles generated by industrial activities (Figure 5(B)). However, as leaded gasoline was phased out in 2000 in Poland [94] the contribution of road traffic may be considered as not significant to the Pb levels in the ŚNP lichens ( Figure 5(B); purple square), which is confirmed by the observed lower 208 Pb/ 206 Pb and higher 206 Pb/ 207 Pb ratios. It follows that we can hypothesise that Pb in our lichen samples is controlled by a binary mixing between terrigenous and industrial Pb. We can then model the mixing relationship for these two endmembers using isotope mass balance equations, and ultimately estimate their corresponding contributions for Figure 4. Relationship between the reciprocal of the organic sulfur concentration (%S org ) and its sulfur isotope composition (δ 34 S org ) in the lichen thalli in our study area. The δ 34 S ranges of regional sources that can potentially contribute are also reported. each sample. The equation of mixing for a k component system for isotope ratios I can be expressed in the following way [95]: with k i=1 z i = 1 and 0 ≤ z i ≤1, and where I is the isotope ratio in the mixture, I i the isotope ratio in each end-member component, z i the mass fraction of component 1, 2, 3 … and n i the relative enrichment of element I in component 1, 2, 3 … relative to component k (i.e. n i = c i /c k ; where c is the concentration of element I). We considered the following characteristics for the two end-members: (i) terrigenous Pb; 206 Pb/ 207 Pb = 1.197 and 208 Pb/ 206 Pb = 2.066 [96], (ii) industrial Pb; 206 Pb/ 207 Pb = 1.155 and 208 Pb/ 206 Pb = 2.112 [96]. The respective contributions calculated for terrigenous and industrial Pb, for each lichen sample, are reported in Table S6.
Results indicate that for both sampling seasons industrial Pb mostly controls (≥60% on average) the budget of this metal in the lichen samples. We are also observing a slightly higher average contribution of industrial Pb during the summer season (63.8 ± 2.0% compared to 60.6 ± 3.2% in winter; Table S6), which is unexpected and requires further studies. Moreover, similarly to the study of Pb concentrations, the highest calculated industrial contributions corresponded to lichen samples collected on peaks. It comforts our previous conclusions that long-range transport of contaminants, probably originating from the Upper Silesian Industrial Region (e.g. smelters located south and west of the study area (Figure 1)), plays a major role in the local air pollution.

Multi-isotope approach
The study site is surrounded by potential sources of air contamination (Figure 1) that include a steel plant, a zinc smelter, sulfur chemical plants, and heat and power stations. The following sections discuss these different emission sources and their typical multiisotope characteristics. Figure 6 reports the coupled S, N and Pb isotope compositions/ratios measured in our lichen samples in an attempt to evaluate the added value of a multi-isotope approach compared to a single one. We separated the summer (black diamond) and winter (white diamond) seasons and focused on the following sources of pollution: A-industrial activities, B-home-heating, and C-road traffic. For each scenario (A, B and C) we selected at least three sampling locations that were expected to be mainly impacted by the corresponding source of pollution.
Results from samples 11, 12, 15, collected on peaks ( Figure 6(A)), suggest that the coupled use of S and Pb isotope systematics brings the most constraints when trying to identify the impact of industrial emissions. Generally, lower 206/207 Pb ratios in the bioindicators indicate a higher contribution of industrial emissions [96], and in particular of coal and galena smelting [97,98]. For δ 34 S, an enrichment in 34 S in winter samples can be observed. While the range of δ 34 S values for summer is wider, for winter it is narrow. This indicates a dominant source of pollution in winter, which is isotopically identified as the industrial activities that are generating SO 4 2-, whether through precipitations [99,100] or/and by the oxidation of industrial SO 2 [101,102]. Figure 6(A) shows that N isotope compositions are very similar but, again, the range is narrower in winter than in summer, indicating additional inputs of other sources of N for the later.
Samples 3, 4, 6, and 7, collected in the vicinity of inhabited areas are most likely impacted by emissions from home heating (in winter) and fertilisers and manure spreading (in summer). In this case, Figure 6(B) suggests that combining S and N isotope systematics helps identifying this type of emission sources. As previously discussed, a decrease in the 206/207 Pb ratio indicates a higher contribution of coal combustion [97,98], whereas more radiogenic ratios may be attributed to inputs of resuspended soil particles [36]. However, here due to too small differences between samples the differentiation was impossible. S isotopes are strongly modified by anthropogenically generated sulfur which can be seen in winter samples. Similarly to the industrial scenario, δ 34 S are controlled by S ions from precipitations in summer and by the oxidation of local SO 2 from home heating in winter. In the case of δ 15 N a change in the isotopic compositions between seasons is also noted, with an enrichment in 15 N in winter that results of higher atmospheric NO x concentrations, generated by coal/wood burning during this season. Figure 6(C) focuses on the eventual influence of emissions from road traffic (diesel and gasoline), using samples from locations 9, 10, 19 that were collected close to roads (Figure 1). It shows that combining the N and S isotope systematics creates a reliable discriminating tool to evaluate the contribution of road traffic to air contamination. As lead has been phased out in Polish gasoline since 2000, seasonal differences in the Pb ratios cannot be explained by a change in the tourist attendance. Also, the ranges of Pb ratios are very narrow and similar between the two seasons. The N isotopes clearly discriminate the sampling seasons and the increase in the summer δ 15 N can be attributed to higher NO x inputs from road traffic during that period. For S isotopes, we are observing a 34 S depletion during the summer while in winter the δ 34 S values are higher. This can be explained by additional influence of other fossil burning sources that enrich lichens in 34 S during winter.

Principal component analysis
Principal component analysis (PCA) is a statistical technique commonly used for identifying patterns in data of high dimension that we implemented here to evaluate the parameters controlling the observed inter-seasonal variability. The difficulty resided in investigating datasets obtained from two distinct media, i.e. bioindicators and the atmosphere. Ultimately, we elected to include the following parameters in the matrix: contaminant's atmospheric concentration, elevation, carbon, nitrogen and sulfur concentrations in the lichens and their corresponding isotopic compositions (δ 13 C, δ 15 N and δ 34 S). All these parameters were also discriminated according to their sampling season: winter vs. summer. When data were missing, the corresponding average value was used instead. Results of the PCA indicate that 80% in winter and 82% in summer of the variability observed in the samples were controlled by five factors (Table 1). The respective remaining 20 and 18% constituted a random noise that is not interpretable using this technique [53,55].
All year long, atmospheric SO 2 concentrations were positively correlated with the elevation, which again confirms that this contaminant is mostly driven by long-range transport (Factor 3 for winter season and Factor 1 for summer season). The highest scores for these two factors corresponded to samples collected on peaks (11, 12, 15; Figure 7) as well as point 13, which is located on the windward side, more exposed to the pollution from long-range transport. Atmospheric CO 2 concentrations and their isotopic compositions (δ 13 C CO2 ) were also correlated (Factor 1 in winter and Factor 2 in summer) but were distinct from the correlation existing between %C org in lichens and Figure 7. Factor scores for lichen samples collected in the Świętokrzyski National Park. PCA factors 1 (summer) and 3 (winter) were selected as they represent the variability associated to the long-range transport of air contaminants.
its δ 13 C org (Factor 3 in summer but no correlation in winter). This confirms our conclusion that δ 13 C org is not controlled by anthropogenic CO 2 emissions but rather by natural sources. The negative correlation between these two parameters during summer may result of an unidentified process controlling the organic carbon in the lichens during this period but not in winter. For nitrogen, we observed a correlation between atmospheric NO 2 and the δ 15 N org in lichens (Factor 4), but only during summer. That hints at a better bioavailability of a 15 N-enriched soil nitrogen during that period. For all seasons, sulfur and nitrogen results showed there is no correlation between the concentrations of these elements and their corresponding isotopic compositions.

Conclusions
We provide here the first report of the multi-isotope characteristics (C, N, S and Pb) of bioindicators collected in the ŚNP in Poland and helped better constrain the compositions and dynamics of regional atmospheric contamination. Results lead to the following conclusions: . The analysis of δ 15 N indicated that coal burning (in winter) and local road traffic with agricultural activities (in summer) are the dominating factors. In the case of δ 34 S, especially at sampling points located on peaks, long-range transport was dominant. This conclusion was confirmed by the PCA analysis. The radiogenicity of the Pb isotope ratios in the studied bioindicators also demonstrated that this area may be impacted by atmospheric emissions of the industrial activities from the Upper Silesian Industrial Region, in agreement with the local topography and dominating winds. Additionally, the presence of the mountain ranges in the centre of the ŚNP acts as a natural barrier and causes that atmospheric emissions from the Upper Silesian are mostly contained within its southern part. On the other hand, carbon isotope compositions were mostly controlled by natural factors that are overprinting the δ 13 C of any atmospheric anthropogenic input. . The seasonality in isotopic compositions of lichens is notable in a few cases, which is due to both the shift in dominating sources between summer and winter and different/preferentially assimilation of the pollutants by lichens. Therefore, the increase of gaseous pollutants (SO 2 , NO 2 ) in winter does not always translate into a visible change in their corresponding δ 15 N and δ 34 S. For Pb ratios differences were not detectable due to a too narrow range of values. . Finally, implementing a multi-isotope approach helps (i) increase the source discrimination power, (ii) identify the main emission sources and/or existing atmospheric secondary processes, and (iii) specifically target the impact of a given source of air contamination by coupling two or three of these isotope systematics. This multiproxy approach can be up-scaled and implemented elsewhere.

Disclosure statement
No potential conflict of interest was reported by the author(s).