In Poland, a national park is a protected area because of its advantages; mainly its outstanding natural value for environmental, scientific, social, cultural, and educational reasons. The national parks are one of the Poland’s strategic natural resources. Therefore, monitoring of the contamination degree of national parks is very important. Industrial and agricultural activities in national parks are prohibited, but airborne chemicals can travel long distances from their sources and therefore can affect ecosystems over broad spatial scales and at locations far from the emission sources (Samecka-Cymerman et al. 2012). Trace elements are emitted by various combustion sources, which can undergo chemical reactions in the atmosphere and fall to the earth as wet deposition (rain, snow), occult deposition (cloud, fog), or dry deposition (dry particles, gas) (Radaman 2004; Schröder et al. 2010). After atmospheric deposition, toxic metals can accumulate in the environmental compartments and tissues of living organisms, increasing the risk of harmful effects (Livett 1988; Seinfeld and Pandis 2006). National parks, which are considered to be the most natural and least contaminated areas, also are affected (Grodzińska 1978).

The problem of pollution by trace elements is still unsolved, and in some cases growing. Therefore, there is a need to continue deposition analysis of these pollutants in the environment (Kosior et al. 2010). Today, moss biomonitoring is a part of pollution monitoring programmes in most European countries as it gives evidence of the anthropogenic impact in urban areas due to vehicular traffic and fossil fuel combustion. It also identifies other sources of heavy metal pollution, such as ore exploitation, agricultural activities, etc. (Aničić 2006; Barandovski et al. 2008; Frontasyeva et al. 2016). Results from moss surveys allow examination of both spatial and temporal trends in metal concentration/deposition, and the identification of areas exposed to the high deposition of metals from long-range transport and local sources. Mosses have a high interception potential for particulate material and precipitation and, therefore, contain significantly higher trace element concentrations than vascular plants and fungi (Ugur et al. 2003; Skuterud et al. 2005). Due to their continuous accumulation of elements, mosses offer information about the sources of pollution long after the pollution episode itself took place (Golubev et al. 2005).

We chose the terrestrial carpet-forming, ectohydric moss Pleurozium schreberi for this study. Ectohydric mosses obtain most trace elements and nutrients directly from the atmosphere (wet and dry deposition) with little uptake from the substrate and therefore are suitable for monitoring airborne pollutants. This species has been demonstrated to be a suitable bioindicator of inorganic substances originating from the atmosphere and has been widely used to map and monitor airborne pollution in European countries (Halleraker et al. 1998). Due to the lack of root systems, mosses depend largely on atmospheric depositions for their nutrient supply and have a great capacity to retain many elements. Mosses lack a cuticle, or it is very thin, due to which their leaves are highly permeable to ions of trace elements, and they usually have a high surface-to-volume ratio (5–10 times higher than vascular plants) (Tyler 1990). Element concentrations as measured in moss samples are supposed to represent the accumulated load of the past 2–3 years (Pakarinen and Rinne 1979). Due to this accumulation, element concentrations measured in moss are much higher than in other sample materials, such as precipitation, dust, or other plants, and thus are easier to measure (Asakawa et al. 2013).

The main objective of the present study was to test the statistical approach of the t quantile to compare elements bioaccumulated in P. schreberi collected from Polish national parks (PNP) and to compare the results between lowland and highland/mountain parks. Our study also examined the view that special care should be taken in regional comparisons of bryophyte concentrations of trace elements. A minor objective is to compare the results with earlier data from Poland. Novel statistical analysis techniques have been used to evaluate the results.

Materials and Methods

Sampling Design

Samples of P. schreberi moss were collected from four lowland parks: Bory Tucholskie NP, Wielkopolski NP, Kampinoski NP, Roztoczański NP, and from five highland/mountain parks: Karkonoski NP, Gór Stołowych NP, Ojcowski NP, Babiogórski NP, and Świętokrzyski NP (Fig. 1). In the figures, these names are abbreviated to BT, Wi, Ka, Ro, Ks, GS, Oj, Bg, and Sw, respectively. Mosses were collected from an area measuring about 50 m × 50 m and at least 5 m away from the canopy of the trees so as not to be directly exposed to through fall precipitation. Within the central part of each of the selected site, three sub-squares of 2 × 2 m, covered with P. schreberi, were selected randomly for the collection of this native moss (only the green parts of moss were taken). Each plant sample consists of a mixture of three subsamples. The design of the moss samples sites was performed according to the protocol adopted within “Heavy Metals in European Mosses: 2000/2001 Survey” (UNECE ICP Vegetation 2011). Sampling and handling was performed using polyethylene gloves and bags. After selection and removing foreign bodies, the moss samples were transported to the laboratory. In the laboratory, the samples were carefully cleaned of all dead material and attached litter.

Fig. 1
figure 1

Localisation of the studied national parks. Samples of P. schreberi moss were collected from four lowland parks: Bory Tucholskie [BT], Wielkopolski [Wi], Kampinoski [Ka], Roztoczański [Ro] and from five upland and mountainous parks: Karkonoski [Ks], Gór Stołowych [GS], Ojcowski [Oj] Babiogórski [Bg] and Świętokrzyski [Sw]

Full size image

Laboratory Analyses

Instrumental Neutron Activation Analysis

The concentration of elements in the P. schreberi moss samples was determined by a multielement instrumental neutron activation analysis (NAA) at the IBR-2 reactor, FLNP JINR, and Dubna. The characteristics of neutron flux density in the two irradiation channels equipped with the pneumatic system and the registration of gamma spectra can be found elsewhere (Frontasyeva 2011). To determine the concentration of elements with short lived isotopes (Al, Ca, Cl, I, Mg, Mn, and V), samples were irradiated for 3 min and measured for 20 min. To determine the concentration of elements with long lived isotopes (Na, Sc, Cr, Fe, Co, Ni, Zn, As, Se, Rb, Sr, Zr, Mo, Sb, Cs, Ba, La, Ce, Sm, Eu, Tb, Hf, Ta, W, Th, and U) (Table 1), the cadmium-screened channel 1 was used. Samples were irradiated for 4 days, repacked, and measured twice using HP germanium detectors after 4 and 20 days of decay, respectively. The NAA data processing and the determination of element concentrations were performed using software developed in FLNP JINR (Dmitriev and Yu 2013).

Table 1 Radionuclides used for neutron activation analysis and their γ-energies
Full size table

Quality Control for Instrumental Neutron Activation Analysis

To evaluate the precision and accuracy of the results, the standard reference materials (SRMs) [1572 (citrus leaves), 2710 (Montana soil), 1632c (trace elements in coal), and 1547 (pPeach leaves)] were simultaneously irradiated with the analysed samples. The results obtained for the SRMs, as well as their certified values, are shown in Table A—Supplementary Information. The measured concentrations and certified values are in good agreement, with the overall precision for measurement between 1 and 10% for most of the elements, except Mo and Sr, which were found to have a higher uncertainty.

Atomic Absorption Spectrophotometry Analysis

The total content of Cd, Cu, and Pb was determined by the ASS method. The unwashed samples were air-dried at 30 °C. To conduct the determination of the elements, the samples were digested with 3 ml of nitric acid (ultra-pure, 65%) and 2 ml of perchloric acid (ultra-pure, 70%) in a CEM Mars 5 microwave oven (Matusiewicz 2003). The plant digests were analysed for Cd, Cu, and Pb using an ETAAS with a Graphite Furnace GF3000, both Atomic Absorption Spectrophotometry AVANTA PM GBC Scientific Equipment (Lajunen and Peramaki 2004). All elements are determined against standards (Atomic Absorption Standard Solution from Sigma Chemical Co.) and blanks containing the same matrix as the samples and subjected to the same procedure. All results for the plants are calculated on a dry weight basis. The accuracy of the methods applied for the determination of elements after microwave-assisted digestion of plants and soil samples is checked by analysing certified reference materials. We used bush branches and leaves from DC73348 LGC standards and RTH 907 Dutch Anthropogenic Soil (Wageningen Evaluating Programmes for Analytical Laboratories, WEPAL) as certified reference materials (Table B—Supplementary Information). The coefficient of variance (CV) is calculated for the determined concentrations of elements in the reference materials.

Statistical Analysis

Because concentration represents the amount of a selected substance in the total amount of all components, its analysis requires application of appropriate mathematical methods (Aitchison 2003; Pawlowsky-Glahn and Buccianti 2011). The sample space of our measurement results is a set of compositions. Composition of a system composed of D components can be calculated by the closure operation C, defined as follows:

$$C({\mathbf{z}}) = \left[ {\frac{{k_{s} {\mathbf{z}}_{1} }}{{\sum\nolimits_{j = 1}^{D} {{\mathbf{z}}_{j} } }},\frac{{k_{s} {\mathbf{z}}_{2} }}{{\sum\nolimits_{j = 1}^{D} {{\mathbf{z}}_{j} } }}, \ldots ,\frac{{k_{s} {\mathbf{z}}_{D} }}{{\sum\nolimits_{j = 1}^{D} {{\mathbf{z}}_{j} } }} = {\mathbf{x}}} \right]$$
(1)

where z is the vector of component amounts (mass or count, for example), k depends on the units of measurement and x is the vector of the compositions. For mg/kg, the k value is 106. The subcomposition can be obtained by applying the closure operation to a vector of S parts taken from z (S < D).

To describe the difference between two compositions, the scalar measure of distance was introduced. The Aitchison distance dA between two compositions x[x1,x2,.,xD] and y[y1,y2,.,yD] is described by the formula (Pawlowsky-Glahn and Buccianti 2011):

$$d_{{\text{A}}} ({\mathbf{x}},{\mathbf{y}}) = \sqrt {\sum\limits_{i = 1}^{D} {\left( {\ln \frac{{x_{i} }}{{g({\mathbf{x}})}} - \ln \frac{{y_{i} }}{{g({\mathbf{y}})}}} \right)}^{2} }$$
(2)

An insight into the data structure can be delivered by cluster analysis methods. Clustering comprises unsupervised methods that can be used to organise data into groups based on similarities among the individual data items. Clustering involves dividing data points into clusters so that items in the same class are as similar as possible. Different types of similarity measures may be used to identify classes; in this work, the distance between compositional points in simplex space, defined in Eq. 2, was used. In hard clustering (non-fuzzy), data are divided into crisp clusters, where each data point belongs to exactly one cluster. In fuzzy clustering, the data points can be assigned to more than one cluster. For each of the points, its membership grade is determined to indicate the degree to which the data points belong to the different clusters. Initial results, assuming the number of clusters, were used in fuzzy clustering analysis.

For each sampling point, the contribution of a cluster was calculated. To estimate the contribution of the kth element in the jth cluster structure, the value of the following expression can be considered:

$$FC_{j} = \sum\limits_{i = 1}^{n} {f_{ji} C_{ki} }$$
(3)

A high value of the chemical element concentration and its high contribution in the cluster produces big values of the FCj vector elements. Low concentrations and contributions produce low values of the vector’s elements. A big FCjk vector element assume the significant contribution of the kth element in the jth cluster structure.

The compositional data are not independent of each other; if content of one of the components increases, the others have to decrease. It is the reason why the results of standard statistical analysis of the relationships between raw components or parts are spoiled by spurious effects (Pawlowsky-Glahn and Buccianti 2011; Pawlowsky-Glahn and Egozcue 2006). The particular properties of compositional data preclude the application of standard statistical techniques on such data in its raw form. For example, variances of constrained, compositional data components are not independent. Variance in a component concentration may not be strictly related to variance in the component abundance, because this variability also is affected by variabilities in the abundances of the other components. As the result of the closure operation, variance of a component and its covariances with the other components are bound by the relationship:

$${\text{cov}} ({\mathbf{x}}_{1} ,{\mathbf{x}}_{2} ) + \cdots + {\text{cov}} ({\mathbf{x}}_{1} ,{\mathbf{x}}_{D} ) = - {\text{var}} ({\mathbf{x}}_{1} )$$
(4)

The application of classical data analysis methods for the statistical analysis of compositional data may lead to false or delusive conclusions. This is why in data interpretation statistical methods designed for compositional data elaboration were used. Particularly, the utilisation of a classical correlation estimator in the statistical analysis of compositional data can lead to delusive conclusions (Filzmoser and Karel 2008). An alternative method, appropriate for covariability estimation in compositional data, is based on the variation matrix T. The elements t of this matrix are variances of the concentration pairs’ log-ratios (Aitchison 2003):

$$t_{jl} = {\text{var}} \left( {\ln \frac{{x_{j} }}{{y_{l} }}} \right),j,\,l = 1, \ldots ,D$$
(5)

The variance of the log-ratio is small when changes in both concentrations follow the same trend, i.e., both values simultaneously increase or decrease. The biggest log-ratio variance is observed when concentrations follow opposite trends. An increase in one variable’s values accompanied by a decrease in the values of the second one produces highly differentiated log-ratios. As the result, the biggest t value is observed. To distinguish between these two types of covariability, the terms “positive” (an increase in both concentrations) and “negative” (an increase in one concentration and a decrease in the second) are used. However, interpretation of a high t is restricted. The actual relationship between abundances of components in a system cannot be reflected by the corresponding t value.

The t values between these characteristics for negative or positive covariabilities describe a random or non-linear relationship between the concentrations. The intermediate values of randomly arranged data can be used in the assessment of the actual t classification.

The distribution function of t can be estimated on the basis of the actual data. Data permutation in the first variable, keeping an unchanged sequence in the second one, produces a number of t values. Most of them correspond to a random arrangement of pairs of variables. A fraction of the lowest t values correspond to concentrations positively related with each other. A fraction of the highest t values correspond to negative covariability between concentrations (Ziembik and Dołhańczuk-Śródka 2015). In this work, the t distribution quantile q0.01 was the upper limit for positive covariability, and q0.99 was the lower limit for a negative one. Although, as mentioned earlier, Pearson’s correlation coefficient r is generally not appropriate for covariability assessment among concentrations, its value was calculated to compare conclusions with the ones related to the t estimates.

In computations, the R language (R Development Core Team 2015) and functions from the libraries, “cluster” (Kaufman and Rousseeuw 2005; Maechler et al. 2014) and “compositions” (van den Boogaart et al. 2014) were used.

Results and Discussion

The calculated statistical parameters of the data are presented in Table 2 and Table D—Supplementary Information. The following parameters are shown: (min—minimal value, q1—lower quartile, median, arithmetic mean, q3—upper quartile, max—maximal value, SD—standard deviation, and MADN—normalised median absolute deviation about the median (Maronna et al. 2006). The MADN parameter can be regarded as a robust alternative to the SD. If no outliers are in the vector of the normally distributed data then the MADN and SD values are similar to each other. Comparison of the median and mean values of the concentration indicates skewness in the data distribution. The most significant difference between the mean and median is observed for Sc. The mean concentration of this element was twice as high as the median. For Hf, Co, Fe, Na, V, Cs, Al, Mg, and Cr, the mean to median ratio exceeds 1.5. For Sb, the lowest ratio of 1.0 was obtained. A ratio lower than 1.2 was observed for As, Ca, Br, K, Mo, and Zn. For the SD and MADN, similar differences occurred. The biggest SD to MADN ratio, calculated for Sc, was 6.4. Ratios higher than 3.0 were calculated for Co, Cr, Na, Fe, and V concentrations. The lowest ratios, not higher than 1.2, were found for Mn, Rb, Sb, Mo, U, and Br. Comparison of the statistical parameters between mountain parks and lowland parks is presented in Table 2. Concentrations of most elements (As, Ba, Ca, Cd, Cl, Cr, Cu, Mo, Pb, Sr, and Zn) were higher (Wilcox test, p < 0.05) in mosses collected from highland and mountainous parks than from lowland parks, with the opposite relationship only for Mn. In the first Polish biomonitoring study with mosses (Grodzińska 1978), most of the highest heavy metal concentrations also were found in P. schreberi and Hylocomium splendens species in highland and mountain parks (Ojcowski, Świętokrzyski and Babiogórski), whereas the lowest concentrations were recorded in maritime and lowland parks (Woliński, Słowiński, and Białowieski). In the southern part of the country, there are mountainous regions, where the rates of atmospheric element deposition may vary locally which is related, among others things, to altitude, which strongly influences the amount of precipitation, the dry deposition of gases and particles, and the wet deposition of cloud or fog droplets (Gerdol and Bragazza 2006). The pattern obtained in this study confirms the results from other European areas. In the northern and eastern Alps (Austria), mosses have been collected from transects along altitudinal gradients on five mountain ranges (Zechmeister 1995). The results showed a remarkable increase of Pb, Cd, Zn, and S concentrations as the altitude rises. High levels of precipitation were strongly correlated with heavy metal deposition, and this seems to be the main source of heavy metal fallout at higher altitudes. Larger amounts of windblown, indigenous particles also were considered for several heavy metals (e.g., V) higher levels (Zechmeister 1995). In the study from north Italy (Gerdol and Bragazza 2006), concentrations of anthropogenic pollutants (especially Cd and Pb) peaked at mid altitude (1400–1800 m) where the frequency of cloud cover was the highest. The authors concluded that the deposition of trace elements by cloud water, which may account for a significant fraction of the total deposition of anthropogenic trace elements, is closely related to the cloud cover frequency (Gerdol and Bragazza 2006). On the other hand, in southern Poland there are many sources of pollution from local industry, which may be partly responsible for the higher concentrations of elements in mountainous areas (Grodzińska et al. 1999; Grodzińska and Szarek-Łukaszewska 2001). The Upper Silesian Ecological Hazard Area is within the southern Poland. This region has been recognised as area of ecological disaster (Wcisło et al. 2002). Additionally, this area is affected by heavy traffic and airborne transboundary pollution from Germany and the Czech Republic (Appleton et al. 2000). It also was demonstrated that a group of elements may have crustal origin in areas of upland and highland (Agnan et al. 2014; Aničić et al. 2007), which is confirmed in this research by the convariability method described below. The values in this study were compared with the reference values from Markert et al. (2014). P. schreberi from both groups of parks were characterised by higher levels of Al, Ce, Cr, and Fe (Table 2). Furthermore, mosses from the highland/mountain parks were characterised by higher levels of Cs, Pb, and Sc. This suggests a significant impact of elements regarded to be of crustal origin on the total content of all studied elements.

Table 2 Minimum (min), maximum (max), mean (avg), and standard deviation (SD) of the concentration (mg kg−1) of elements in P. schreberi from national parks divided into two groups: mountainous parks and lowland parks
Full size table

The values were compared to the average concentrations of elements (Cd, Cr, Cu, Fe, Pb, and Zn) for the whole area of Poland. In the case of lowland PNP, we obtained lower values for all elements, with the most significant decreases observed for Cu and especially for Pb (Grodzińska et al. 1999; Table 2). Interestingly, for the mountain parks in this study, we received lower values for Cu and Pb (the most significant decrease being for the latter), the level of Cd was similar for both mountain and lowland parks, whereas for other elements, their concentrations were higher in mountain parks. The concentrations of Cr, Fe, and Pb were lower in comparison with the concentrations for PNP, both mountain and lowland, in 1975 (Grodzińska et al. 1999; Table 2). Again, the most significant decrease was recorded for Pb in both lowland and mountain parks. Since the 1970s, a continuous decrease in the concentrations of heavy metals determined in mosses was observed, which is confirmed by Grodzińska et al. (1999) and by measuring metal emissions by the Central Statistical Office – GUS (GUS 2016) (Table C—Supplementary Information). These authors revealed that the decrease over 20 years (1975–1995) in 12 PNP, which were fairly evenly distributed across Poland, was also evident. Moreover, in the study by Kłos et al. (2015), the analyses of the results of the biomonitoring research carried out during the years 1975–2014 indicate an improvement in the quality of the environment in the studied areas of Poland. The mean concentrations of Ni, Cu, Zn, Cd, and Pb in the samples collected in 2014 in the woodlands of north-eastern Poland were even comparable to the mean concentrations of metals determined in mosses collected in the Svalbard archipelago (N Norway), regarded as relatively clean, during the last three decades. To assess the grouping of trends in the data, clustering methods were applied and the results are presented in a tree diagram (Fig. 2). This is a dendrogram constructed using a divisive algorithm. The branches in the tree are marked with the abbreviations for the names of each national park. A group of branches representing the composition of samples collected in the Babiogórski, Ojcowski, and Góry Stołowe PNP appears in the dendrogram. The composition of the samples from Karkonoski PNP seems to form a cluster, but intrusions of branches representing samples from the other PNPs are observed in its structure. The structure of dendrogram supposes an influence of site localisation on sample composition but it is not critical. Some other, currently uncontrolled, factors significantly influenced the elemental composition of the samples. To look for actually overlapping data structures, the fuzzy clustering method was applied. The algorithm applied in the “fanny” function used in computations produced mostly inconsistent results. Only for the number of clusters = 2 and the membership exponent = 2 were the computation results reasonable. To estimate the influence of an element on the cluster formation, the components of the FC matrix were calculated. A significant contribution of an element in the cluster results in a big component value compared with the components representing other clusters. The components of the FC matrix calculated using our data were similar. The lowest component ratio was 0.8 for the alkali metals Rb and Cs, and the biggest one was 1.4 for Sc and Fe (crustal origin elements mostly). These values assume, at most, a weak trend to cluster formation in the data. Weak grouping in the data enables application of covariability analysis. To assess the linear covariability between pairs of element concentrations, the quantiles of the t parameter distribution were calculated. In Table 3, the results of the computations are shown. Symbols describe the covariability character, i.e., “r” is for random or no linear relationships between concentrations, “p” is for positive covariability and “n” is for negative covariability. For comparison, the commonly used data interpretation Pearson's correlation coefficients r were calculated. The r parameters were tested for a 0 value in the population, at a p-level of 0.02. For 2D normally distributed data composed of variable pairs representing 30 cases, the critical, absolute |rc| value is 0.423. In Table 3, positive correlations bigger than r are in bold and negative correlations less than r are underlined. Linear covariabilites estimated with the t quantiles and the correlation coefficient r often lead to similar, but not identical, conclusions. For example, both the t quantile and r for the Na and Mg pair suppose positive covariability in these metals’ concentrations. But though r for the Al and Na concentrations is big, the t quantile indicates no linear relationship between these variables. In contrast, the poor correlation between Na and K concentrations is associated with the positive covariability supposed by the t quantile. In the assessment of the covariabilities in the data, the t quantile is preferred due to its compositional coherence. For this reason, this parameter is used in statistical inference, and the r parameter influence on the data interpretation is limited.

Fig. 2
figure 2

Clustering methods were applied and the results are presented in the tree diagram. The structures of dendrograms constructed using a divisive algorithm are shown

Full size image
Table 3 Covariabilities in concentration pairs estimated with Pearson's correlation coefficient (above the table diagonal) and the t quantile encoded “r” for random or nonlinear, “p” for positive and “n” for negative covariability
Full size table

Negative covariabilities (increase–decrease) in concentrations were also observed, in most cases for pairs of elements where one of the components was Mn, Rb, or Cd. In this study, manganese in a pair with Zn, Mo, Ba, Pb, and Cd appears to be a potent antagonist of these elements. This supports the thesis by Varela et al. (2015) that the concentrations of the elements can be attributed to causes other than the inherent physicochemical characteristics of the moss. Antagonisms and interactions among the some of these elements (Mo, Zn) and their effects on higher plant growth have been reported by a number of investigators (2011). Mn is the one of elements for which there was no significant correlation between the concentrations in moss and in the bulk deposition (Boquete et al. 2011). The relationship between exposure and concentration is altered by some as yet unknown reason (i.e., the physicochemical characteristics of the emissions, physicochemical and biological processes, and their interactions with environmental factors) (Boquete et al. 2011). Mn is an essential element for moss and, therefore, moss has a high capacity to regulate the levels of this element, unlike those of other contaminants, which may explain the dependence of elements considered to be derived primarily from anthropogenic sources (e.g., Pb and Cd). Although, at high concentrations Mn also may cause toxic effects (El-Joual and Cox 1998). For the alkali metal Rb, we observed negative covariabilities with Cd and Pb. Most probably, Rb can replace Cd and Pb cations in the exchange centres of mosses. For Rb there also were negative covariabilities with Br and Sb.

In most cases, the lowest t value in the complete data (positive covariabilities, increase- increase) was found for pairs of elements released into the environment primarily from anthropogenic (e.g., for pairs of elements regarded as anthropogenic in origin, mostly Pb–Cd, Zn–Cd, As–Sb) and crustal (e.g., for pairs of elements regarded as crustal in origin, mostly Al–V, Sc–Fe, Fe–Co, Th–U) sources. Low values of t assume proportionality of the listed element concentrations, assuming that in moss their constant ratio is required. The statistical methods used in this study can be useful in similar biomonitoring projects, especially in regional comparisons of mosses concentrations of trace elements.

Conclusions

Moss samples were collected in unpolluted areas of national parks. The distribution of element concentrations was differentiated. Some of the distributions were significantly skewed and in the data outliers are expected. Although their elemental composition was not uniform, clear spatial patterns in concentration distributions were not observed. In the linear covariability estimation, the t quantile approach was used. A number of positive covariabilities were observed. This could be a result of the geochemical characteristics of the local environment, including the crust composition to which soil composition is related.

Due to the protection of the national park areas, the influence of pollution sources on the elemental composition of moss is limited. The concentrations of most elements (As, Ba, Ca, Cd, Cl, Cr, Cu, Mo, Pb, Sr, and Zn) were higher in mosses collected from highland and mountainous parks than in lowland parks, with the opposite relationship only for Mn. A comparison of data obtained from Polish national parks in 1970s and 1990s showed a significant decrease in the concentrations of heavy metals.