Atmospheric Heavy Metal Deposition in North Macedonia from 2002 to 2010 Studied by Moss Biomonitoring Technique

: Moss biomonitoring technique was used for a heavy-metal pollution study in Macedonia in the framework of the International Cooperative Program on E ﬀ ects of Air Pollution on Natural Vegetation and Crops (UNECE IPC Vegetation). Moss samples ( n = 72) were collected during the summers of 2002, 2005, and 2010. The contents of 41 elements were determined by neutron activation analysis, atomic absorption spectrometry, and inductively coupled plasma atomic emission spectrometry. Using factor and cluster analyses, three geogenic factors were determined (Factor 1, including Al, As, Co, Cs, Fe, Hf, Na, Rb, Sc, Ta, Th, Ti, U, V, Zr, and rare-earth elements–RE; Factor 4 with Ba, K, and Sr; and Factor 5 with Br and I), one anthropogenic factor (Factor 2, including Cd, Pb, Sb, and Zn), and one geogenic-anthropogenic factor (Factor 3, including Cr and Ni). The highest anthropogenic impact of heavy metal to the air pollution in the country was from the ferronickel smelter near Kavadraci (Ni and Cr), the lead and zinc mines in the vicinity of Makedonska Kamenica, Probištip, and Kriva Palanka in the eastern part of the country (Cd, Pb, and Zn), and the former lead and zinc smelter plant in Veles. Beside the anthropogenic inﬂuences, the lithology and the composition of the soil also play an important role in the distribution of the elements.


Introduction
Elements and compounds that are released into the environment and damage the biosphere (soil, water, air) are treated as pollutants-regardless of whether they come from natural processes on the Earth's surface or human activities [1][2][3]. Atmospheric emissions, due to their widespread dispersion and exposure to living organisms, pose a major risk to human health [4,5]. The representation of pollutants is mostly dependent on the source of emissions, the amount of the released substances, their composition, state of matter, altitude, and weather conditions season, amount of atmospheric sediments, air currents, precipitation, acidity, temperature, humidity; the rate of deposition depends on several factors, such as meteorology wind velocity, relative humidity, and temperature; particle shape and size, and the chemical forms of the elements [1,2,6,7]. The Republic of North Macedonia has three climatic types, due to specific natural and geographical features. The northern parts of the country are characterized by the mildly continental climate, while the area along the valleys of the Vardar River typically has a Mediterranean climate. There is also a mountainous climate, which is dominant in the mountainous region of the country. All three types of climate that spread and partly intertwine through the country have an impact on the properties of the individual seasons between regions. Winters are milder in the eastern areas, and summers are hotter and drier in the western areas [53].
The average annual precipitation varies from 500 to 600 mm in the central and eastern areas along the valley of the Vardar River to 1400 mm in the western mountainous area. The area with the least precipitation is between the Ovče Pole, Tikveš, and the Štip basins. Usually, the wind blows from the north/northwest toward the south/southeast or vice versa. However, the directions of the winds in various parts of Macedonia are heavily dependent on the orographic conditions. The average speed of the winds in the country is 6 m/s [53].
The country is enriched by reserves of several ores. The industrial sector represents 31.9 of Gross Domestic Product (GDP), while the agriculture sector represents 18.9%. The main industrial activities are situated in Skopje (steel and ferroalloys production and steel processing), Radoviš (copper mine and flotation), Veles (abandoned lead-zinc and cadmium smelter), Kavadarci (ferronickel smelter), Tetovo (abandoned ferrochrome smelter and present ferrosilicon smelter), and Kičevo and Bitola (thermoelectric power plants, where the lignite fuel is used). Three Pb-Zn mines and flotation plants are located (Sasa, Toranica, Zletovo) in the eastern part of the country; a previously active As-Tl-Sb mine "Allchar" is in the southern part of Macedonia. These active and abandoned industrial facilities have led to environmental pollution in the nearby and distant areas with different heavy metals, reported by several studies [47][48][49][50] and those studies realized recently [46,[54][55][56][57][58][59].
Macedonia is constituted by several geotectonic units that were formed during different geological periods; these are generally spread out in NNW-SSE to N-S direction [44,59]. Contacts of units are marked by regional faults managed by the Laramide compression (Laramide orogeny phase) [44]. The country lies on four major tectonic units ( Figure 2): The Serbo-Macedonian massif (SMM), the Vardar zone (VZ), the Pelagonian massif (PM), and the West-Macedonian zone (WMZ) The average annual precipitation varies from 500 to 600 mm in the central and eastern areas along the valley of the Vardar River to 1400 mm in the western mountainous area. The area with the least precipitation is between the Ovče Pole, Tikveš, and the Štip basins. Usually, the wind blows from the north/northwest toward the south/southeast or vice versa. However, the directions of the winds in various parts of Macedonia are heavily dependent on the orographic conditions. The average speed of the winds in the country is 6 m/s [53].
The country is enriched by reserves of several ores. The industrial sector represents 31.9 of Gross Domestic Product (GDP), while the agriculture sector represents 18.9%. The main industrial activities are situated in Skopje (steel and ferroalloys production and steel processing), Radoviš (copper mine and flotation), Veles (abandoned lead-zinc and cadmium smelter), Kavadarci (ferronickel smelter), Tetovo (abandoned ferrochrome smelter and present ferrosilicon smelter), and Kičevo and Bitola (thermoelectric power plants, where the lignite fuel is used). Three Pb-Zn mines and flotation plants are located (Sasa, Toranica, Zletovo) in the eastern part of the country; a previously active As-Tl-Sb mine "Allchar" is in the southern part of Macedonia. These active and abandoned industrial facilities have led to environmental pollution in the nearby and distant areas with different heavy metals, reported by several studies [47][48][49][50] and those studies realized recently [46,[54][55][56][57][58][59].
Macedonia is constituted by several geotectonic units that were formed during different geological periods; these are generally spread out in NNW-SSE to N-S direction [44,59]. Contacts of units are marked by regional faults managed by the Laramide compression (Laramide orogeny phase) [44]. The country lies on four major tectonic units ( Figure 2): The Serbo-Macedonian massif (SMM), the Vardar zone (VZ), the Pelagonian massif (PM), and the West-Macedonian zone (WMZ) [60,61]. The SMM and PM form the oldest complexes in which highly metamorphic Proterozoic rocks have developed; the area is characterized by widespread Riphean-Cambrian, Paleozoic, Mesozoic, and Cenozoic rocks [44]. The WMZ is lithologically composed of low-grade metamorphic rocks, and anchi-metamorphic Paleozoic rocks and migmatites, Triassic and Jurassic sediments, and migmatites, and Tertiary sediments [44]. The VZ is a large and important lineament structure of the Balkan Peninsula [61]; it represents an area where a continuous thinning of the continental crust had occurred, and where, during the Lower Jurassic, it completely transited into the rift zone [44].
Atmosphere 2020, 11, x FOR PEER REVIEW 4 of 24 [60,61]. The SMM and PM form the oldest complexes in which highly metamorphic Proterozoic rocks have developed; the area is characterized by widespread Riphean-Cambrian, Paleozoic, Mesozoic, and Cenozoic rocks [44]. The WMZ is lithologically composed of low-grade metamorphic rocks, and anchi-metamorphic Paleozoic rocks and migmatites, Triassic and Jurassic sediments, and migmatites, and Tertiary sediments [44]. The VZ is a large and important lineament structure of the Balkan Peninsula [61]; it represents an area where a continuous thinning of the continental crust had occurred, and where, during the Lower Jurassic, it completely transited into the rift zone [44].

Sampling and Sample Preparation
In this study, datasets from three different moss surveys were used for comparison among them. The first moss survey aiming to monitor heavy metal air pollution was carried out in 2002 in the entire territory of Macedonia. For this reason, samples from 72 sites of three moss species Hypnum cupressiforme, Homalothecium lutescens, and Homalothecium sericeum were collected from September to October [1,2]. The contents of 41 elements (Na, Mg, Al, Cl, K, Ca, Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, As, Se, Br, Rb, Sr, Zr, Mo, Ag, Cd, In, Sb, I, Cs, Ba, La, Ce, Sm, Eu, Tb, Hf, Ta, W, Au, Th, and U) in each sample were analyzed either by neutron activation analysis or by atomic absorption spectrometry. The next survey was conducted from August to September 2005; samples of Hypnum cupressiforme and Homalothecium lutescens were taken at the same 72 sites. Using neutron activation analysis and atomic absorption spectrometry, the content of 38 elements (Al, As, Ba, Br, Ca, Cd, Ce, Co, Cr, Cs, Cu, Dy, Eu, Fe, Hf, Hg, K, La, Mg, Mn, Mo, Na, Ni, Rb, Sb, Sc, Se, Sm, Sr, Ta, Tb, Th, U, V, W, Yb, Zn, and Zr) were determined in each sample [32]. The third moss survey was conducted in August and September 2010 at the same 72 locations, where the same moss species were collected as in the moss survey carried out in 2005. In addition to the techniques used in previous surveys, the atomic emission spectrometry with inductively coupled plasma (ICP-AES)

Sampling and Sample Preparation
In this study, datasets from three different moss surveys were used for comparison among them. The first moss survey aiming to monitor heavy metal air pollution was carried out in 2002 in the entire territory of Macedonia. For this reason, samples from 72 sites of three moss species Hypnum cupressiforme, Homalothecium lutescens, and Homalothecium sericeum were collected from September to October [1,2]. The contents of 41 elements (Na, Mg, Al, Cl, K, Ca, Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, As, Se, Br, Rb, Sr, Zr, Mo, Ag, Cd, In, Sb, I, Cs, Ba, La, Ce, Sm, Eu, Tb, Hf, Ta, W, Au, Th, and U) in each sample were analyzed either by neutron activation analysis or by atomic absorption spectrometry. The next survey was conducted from August to September 2005; samples of Hypnum cupressiforme and Homalothecium lutescens were taken at the same 72 sites. Using neutron activation analysis and atomic absorption spectrometry, the content of 38 elements (Al, As, Ba, Br, Ca, Cd, Ce, Co, Cr, Cs, Cu, Dy, Eu, Fe, Hf, Hg, K, La, Mg, Mn, Mo, Na, Ni, Rb, Sb, Sc, Se, Sm, Sr, Ta, Tb, Th, U, V, W, Yb, Zn, and Zr) were determined in each sample [32]. The third moss survey was conducted in August and September 2010 at the same 72 locations, where the same moss species were collected as in the moss survey carried out in 2005. In addition to the techniques used in previous surveys, the atomic emission spectrometry with inductively coupled plasma (ICP-AES) was also used to determine 41 elements (Al, As, Au, Ba, Br, Ca, Cd, Ce, Cl, Co, Cr, Cs, Cu, Dy, Eu, Fe, Hf, Hg, I, K, La, Mg, Mn, Mo, Na, Ni, Pb, Rb, Sb, Sc, Se, Sm, Sr, Ta, Tb, Th, Ti, U, V, W, Zn, and Zr) in the samples [22,33].
To enable a comparison between individual surveys, the sampling locations mostly overlap (Figure 3), and the same biotope conditions were used. The sampling in 2002 followed the principles of European moss surveys [14,62]. In 2005, sampling was performed according to the guidelines of the ICP Vegetation Program-monitoring manual for the 2005/2006 survey, and the procedure used in the previous European moss surveys [63]. In 2010, the sampling followed the principles of the Convention on Long-Range Transboundary Ari Pollution (CLRTAP) and the International Cooperative Program on Effects of Air Pollution on Natural Vegetation and Crops monitoring manual for 2010/2011 survey and the procedure used in the previous European moss surveys [64].
To enable a comparison between individual surveys, the sampling locations mostly overlap (Figure 3), and the same biotope conditions were used. The sampling in 2002 followed the principles of European moss surveys [14,62]. In 2005, sampling was performed according to the guidelines of the ICP Vegetation Program-monitoring manual for the 2005/2006 survey, and the procedure used in the previous European moss surveys [63]. In 2010, the sampling followed the principles of the Convention on Long-Range Transboundary Ari Pollution (CLRTAP) and the International Cooperative Program on Effects of Air Pollution on Natural Vegetation and Crops monitoring manual for 2010/2011 survey and the procedure used in the previous European moss surveys [64]. In each sampling year, the samples were collected at least 300 m from a major road and populated area, and at least 100 m from small roads and single houses. Each sample was a mix of five to ten sub-samples taken from 50 m × 50 m. Sampling and sample handling was performed using separate polyethylene gloves for each sample, to prevent any contamination of the samples. Samples were properly marked and stored in paper bags. In cases where more than one moss species was collected at the same sampling site, the interspecies comparison showed no essential differences within error estimates.
Samples were brought to the laboratory where they were cleaned from extraneous material (litter and dead leaves) and dried to a constant weight for 48 h at 30-40 °C. Green-brown parts of the moss shoots represented the last three years of the moss growth and were selected to be subjected to the analysis of pollutant content. For neutron activation analysis (NAA), atomic absorption spectrometry (AAS), and atomic emission spectrometry with inductively coupled plasma (ICP-AES), the moss samples were prepared according to the description from previous studies, where quality control data are also given [2,22,32]. In each sampling year, the samples were collected at least 300 m from a major road and populated area, and at least 100 m from small roads and single houses. Each sample was a mix of five to ten sub-samples taken from 50 m × 50 m. Sampling and sample handling was performed using separate polyethylene gloves for each sample, to prevent any contamination of the samples. Samples were properly marked and stored in paper bags. In cases where more than one moss species was collected at the same sampling site, the interspecies comparison showed no essential differences within error estimates.
Samples were brought to the laboratory where they were cleaned from extraneous material (litter and dead leaves) and dried to a constant weight for 48 h at 30-40 • C. Green-brown parts of the moss shoots represented the last three years of the moss growth and were selected to be subjected to the analysis of pollutant content. For neutron activation analysis (NAA), atomic absorption spectrometry (AAS), and atomic emission spectrometry with inductively coupled plasma (ICP-AES), the moss samples were prepared according to the description from previous studies, where quality control data are also given [2,22,32].

Statistical Methods
Data on chemical analyses obtained from previous studies [1,2,22,32,33] are used in the proposed report. To compare the results between different datasets, and simultaneously perform parametric and non-parametric statistical methods, the statistical software Statistica 13 (StatSoft, Inc., Tulsa, OK, USA) was used [65,66].
Since the datasets were non-normally distributed, a Box-Cox transformation was performed [67]. Box-Cox transformed values were used later in the analysis of variance (ANOVA). This analysis represents an important tool in assessing the differences in selected elements according to the sampling campaign. With the assist of bivariate statistics, the correlation of content of chemical elements between individual moss samples according to the sampling campaign was initially established, which was then substantiated by statistical tests (t-test and F-test). The use of t-test and F-test allows the assessment of the statistical significance of differences between the accumulation values measured in different exposure periods. The enrichment ratio (ER) was used to assess the selected elements measured in various sampling campaigns. The ER is defined as the ratio of a grade of an element in some areas according to the natural abundance of the metal [68].
To investigate the degree of association of the chemical elements in the moss samples, the Pearson correlation coefficient I was used [69]. Using the correlation coefficient, the correlation degree (linear dependence) between two random variables or sets of random variables was calculated. A good association between the elements was assumed for r between 0.5 and 0.7, while a strong association was determined if r values ranged between 0.7 and 1.0. The elements with small values of the correlation coefficient were omitted from the matrices.
Based on the matrix of correlation coefficients, multivariate cluster, and R-mode factor analysis, the associations of chemical elements were revealed [70][71][72]. In this case, the values in a range of 0.5 to 0.7 were a good association, while the values in the range of 0.7 to 1.0 were a strong association. The multivariate statistical cluster and factor analyses were performed on 27 selected elements (Al, As, Co, Cs, Fe, Hf, Na, Rb, Sc, Ta, Th, Ti, U, V, Zr, Cd, Pb, Sb, Zn, Cr, Ni, Ba, K, Sr, Br, I, and rare-earth elements (RE)). Some elements were excluded from the analysis because they did not show a sensible connection with other chemical elements, or did not satisfy the criteria of dimension variables based on several observations. For factor analysis purposes, the variables were standardized to a zero mean [70,71,73]. Subsequently, based on characteristics of the 27 individual elements, five synthetic variables were designed (Factor 1-F1 to Factor 5-F5), which account for 81.8% of the total variability of the treated elements. The variables with factor loadings higher than 0.6 were assumed to contribute significantly to a given factor.
For the construction of maps showing the spatial distribution of factor scores and maps displaying the distribution of heavy metals in moss samples, universal kriging with the linear variogram interpolation method was applied [70]. For interpolation, the basic grid cell size was 1 × 1 km. Seven classes (0-10, 10-25, 25-40, 40-60, 60-75, 75-90, and 90-100) of the percentile values of the distribution of interpolation values were defined. The visualization (mapping) of data was performed using several software packages: Statistica 13 (StatSoft, Inc., Tuls, OK, USA), QGIS (#), and Surfer 17 (Golden Software, Inc., Golden, CO, USA).

Results and Discussion
The results of descriptive statistics to all assessed elements (aluminum, arsenic, barium, bromine, cadmium, cerium, cobalt, chromium, cesium, europium, iron, hafnium, iodine, potassium, lanthanum, sodium, nickel, lead, rubidium, antimony, scandium, samarium, strontium, tantalum, terbium, thorium, titanium, uranium, vanadium, zinc, and zirconium) in the 72 samples are shown in Table 1. The mean, median, minimum, and maximum contents were calculated for all 31 chemical elements (Al, As, Ba, Br, Cd, Ce, Co, Cr, Cs, Eu, Fe, Hf, I, K, La, Na, Ni, Pb, Rb, Sb, Sc, Sm, Sr, Ta, Tb, Th, Ti, U, V, Zn, and Zr) for the three sets of samples ( Table 1). The results of the ANOVA and ER with respect to the sampling campaign are given in Table 2.  In Table 1, the mean values, median values, and the interval ranges for the contents of all 31 selected elements used for comparison between moss surveys measurements are introduced. Generally, the medians of the elements decreased according to the sampling campaign. However, there are some exceptions. The median content of Ni content from 2005 (5.8 mg/kg) and 2010 (4.3 mg/kg) was higher than in 2002 (2.5 mg/kg), which is related to reactivation and increasing the production capacity of the ferronickel smelter near Kavadarci in 2004 [33].
A The presence of iodine and bromine as an association of these two elements (halogens) is primarily due to marine influence. Along the valley, the Vardar River, the Jugo wind blows from the northern parts of Africa across the Aegean Sea and brings with itself particles enriched with these two elements, which could be subsequently trapped in mosses [74]. The titanium also does not follow a declining trend over the years. The values of Ti from 2005 are higher than in 2002, which can be related to the increased mining activities in the copper flotation plant near Radoviš [22].
Al, As, Ba, Ce, Co, Cs, Eu, Fe, Hf, La, Na, Rb, Sc, Sm, Ta, Tb, Th, U, V, and Zr represent the typical crustal composition. However, it cannot be overlooked the significant fall of some median values of elements according to the sampling campaign. One possible explanation in this case study could be related to the intensity of precipitation and leaching of trapped particles from the mosses, reducing the values by as much as 20%, according to previous studies [75]. In eight regions in Macedonia, the precipitation increased between 2005 and 2010, which further confirms the aforementioned study [76,77].
The differences in the measurements of individual elements, according to the sampling campaign, were supported by the ANOVA analysis. The t-test was used to compare the values of the elements obtained in two different sampling campaigns (Table 2). There were no significant differences between the contents of particular elements measured in 2002 and 2005, which means that all mining and metallurgy activities continuously operated within the same capacity during that period; there were also no natural changes. There were significant differences between the contents of elements measured in 2002 vs. 2010 and the contents of elements obtained in 2005 vs. 2010, respectively, related to either natural or anthropogenic influences. This was then also confirmed by the F-test, which showed that there are significant differences in element values in all three moss surveys.
The enrichment ratio (ER) for all selected 31 elements is given in Table 2. The highest enrichment ratio is found for Ni (2.34) in moss samples from 2002 compared with those from 2005 (Table 2). High values of Ni were a consequence of the increased workload in ferronickel smelter in the vicinity of Kavadarci in 2004. The highest enrichment ratio was found for Br (1.96), which is a marine halogen. The Aegean Sea is less than 60 km from the southern state border, and the wind blowing from the northern parts of Africa across the Aegean Sea brings halogen-enriched particles. Over the years, the wind may be strengthened, increasing the content of Br in mosses [78].
For a better visibility, the matrix of correlations coefficients was divided into two matrices. Table A1 is a matrix of correlation coefficient of 16 selected naturally distributed elements; Table A2 is a matrix of correlation coefficients of 11 selected elements. In both matrices, the Box-Cox transformed values were used.
Three geogenic, one anthropogenic association, and one mixed (geogenic-anthropogenic) geochemical associations were established as a result of a visual inspection of the similarities of the spatial distribution of element patterns, a comparison of basic statistical parameters, the correlation coefficient matrices (Tables A1 and A2), and the results of multivariate analyses (cluster and factor analysis). Variables or samples with similar behavior were combined into the cluster according to the Pearson correlation coefficient. The results are presented in a dendrogram, which was performed to show the results of the hierarchical cluster analysis of 27 selected elements ( Figure 4); and with factor analysis, where the grouping of the elements into five groups was performed ( Table 3). The same classification trend into the factor score or cluster can be seen in both analyses. Four elements were excluded from the clustering because of the absence of similarity linkage. transformed values were used.
Three geogenic, one anthropogenic association, and one mixed (geogenic-anthropogenic) geochemical associations were established as a result of a visual inspection of the similarities of the spatial distribution of element patterns, a comparison of basic statistical parameters, the correlation coefficient matrices (Tables A1 and A2), and the results of multivariate analyses (cluster and factor analysis). Variables or samples with similar behavior were combined into the cluster according to the Pearson correlation coefficient. The results are presented in a dendrogram, which was performed to show the results of the hierarchical cluster analysis of 27 selected elements ( Figure 4); and with factor analysis, where the grouping of the elements into five groups was performed ( Table 3). The same classification trend into the factor score or cluster can be seen in both analyses. Four elements were excluded from the clustering because of the absence of similarity linkage. The matrix of factor loadings is shown in Table 3. Combined, the five factors explain 81.8% of the variability of the threatening elements. Factor 1 (F1) represents the strongest factor, with 47.3% of the total variability. This group of elements represents the geogenic group of elements. Factor 2 The matrix of factor loadings is shown in Table 3. Combined, the five factors explain 81.8% of the variability of the threatening elements. Factor 1 (F1) represents the strongest factor, with 47.3% of the total variability. This group of elements represents the geogenic group of elements. Factor 2 (F2) is the second strongest factor, with 10.4% of the total variability. This factor is associated with elements Pb, Zn, and Cd and represents the anthropogenic elements. The third factor (F3) explains 8.6% of total variability and includes elements of Cr and Ni, whose distribution can be natural or anthropogenic. Factor (F4) represents 8.5% of total variability and includes elements that show natural distribution. Factor 5 (F5) explains 7% of total variability and includes Br and I. Association of Al, As, Co, Cs, Fe, Hf, Na, Rb, Sc, Ta, Th, Ti, U, V, Zr, and RE (Factor 1) represents the typical crustal components that are significantly influenced by the mineral particles that are carried into the atmosphere and later trapped in the moss samples by the wind. Generally, the highest content of these elements was detected in samples collected in the vicinity of Precambrian and Paleozoic shales. Shales decomposed under the influence of various weather conditions, and the minerals of these elements are released into the environment. The annual precipitation is the lowest in the regions of Tikveš, Štip, and Ovče Pole, which leads to erosion of decomposed shales by wind and significantly affects the spread of mineral particles [6]. In other words, the distribution of elements from these groups is independent of urban and industrial areas, but depended on the geological background. The median values of some of these elements have changed over the sampling campaign (Table 1), mostly decreasing, which can also be seen from the spatial distribution of Factor 1 (Figures 5 and 6). This can be supported by the presence of Al in this group. Aluminum compounds are insoluble in water, and the total content of Al in various biological systems originates from dust pollution [79]. Similarly, the presence of iron in the samples originated from dust particles. The highest content of Fe was found in samples collected from the central part of the country, in an area with the smallest annual precipitation and low vegetation, where the land is susceptible to various erosion factors. For the same reasons, the influence of the steel smelter in Skopje on the content of these elements in the samples was not to be observed. insoluble in water, and the total content of Al in various biological systems originates from dust pollution [79]. Similarly, the presence of iron in the samples originated from dust particles. The highest content of Fe was found in samples collected from the central part of the country, in an area with the smallest annual precipitation and low vegetation, where the land is susceptible to various erosion factors. For the same reasons, the influence of the steel smelter in Skopje on the content of these elements in the samples was not to be observed. There was a declining trend in the content of elements grouped into Factor 1, according to the sampling campaign ( Figure 5). In 2002, the highest concentrations of the Factor 1 elements were found in moss samples, collected in the regions of Bitola, Strumica, Radoviš, and Skopje. The high values of these elements (especially uranium) in the samples collected in the Serbo-Macedonian massif and the Pelagonia massif ( Figure 6) can be explained by the existence of uranium ore deposits in these regions [44,59]. The second peak of these elements (mostly related to U and Th) was detected in the samples taken in the southwestern part of the country near Bitola, which can be related to fly ash emissions from a thermoelectric power plant in Bitola, where lignite was used as fuel. It is also related to transboundary pollution from several thermoelectric power plants in northern Greece [2]. The uranium-rich granite deposits in the Strumica region can also affect the presence of high values of uranium [44,59]. During the exploitation of these granites, uranium- There was a declining trend in the content of elements grouped into Factor 1, according to the sampling campaign ( Figure 5). In 2002, the highest concentrations of the Factor 1 elements were found in moss samples, collected in the regions of Bitola, Strumica, Radoviš, and Skopje. The high values of these elements (especially uranium) in the samples collected in the Serbo-Macedonian massif and the Pelagonia massif ( Figure 6) can be explained by the existence of uranium ore deposits in these regions [44,59]. The second peak of these elements (mostly related to U and Th) was detected in the samples taken in the southwestern part of the country near Bitola, which can be related to fly ash emissions from a thermoelectric power plant in Bitola, where lignite was used as fuel. It is also related to transboundary pollution from several thermoelectric power plants in northern Greece [2]. The uranium-rich granite deposits in the Strumica region can also affect the presence of high values of uranium [44,59]. During the exploitation of these granites, uranium-reach dust particles may be released into the environment in considerable amounts, trapped by the moss [2].  The high representation of these elements, especially Al, Sc, Ti, V in the vicinity of Radoviš, can be related to the activities of the copper mine and flotation plant Bučim [27,28,37,80]. High values of these elements (Sb, Fe, and Co) near Tetovo and Kavadraci are related to the previously active ferrochrome smelter in Tetovo and current activity near Kavadarci (ferronickel smelter) [2,79,81].
In 2005, the anomalies were no longer so clear. The anomaly can be seen along the Radika River in Western Macedonia where Neogene clastic sediments predominate [44]. Higher contents of elements from that group were also found in the vicinity of Galičica, due to the Mesozoic carbonate rocks that occur in these areas. Furthermore, areas with high contents of the element group in Factor 1 were found in the area south of Skopje, and in the vicinity of Prilep both areas with Proterozoic metamorphic rocks from the Pelagonian massif [8]. The median values are mostly lower than the median values of the same elements measured in 2002.
The median values for these elements were lowest in 2010 ( Figure 5). According to the State statistical office in 2011 and 2013, the precipitation in eight regions in Macedonia was greater in 2010 than in 2005. However, the central part of the country is characterized by low plant cover and the lowest amount of precipitation, which means these areas are most exposed to erosion, and suspension of soil material [22].
The spatial distribution of Factor 2 (Cd-Pb-Sb-Zn) was not influenced by lithological background and is mainly connected with a lead-zinc smelter (Figure 7). Factor 2 was directly related to the dust from the flotation tailings at the mines and from the soils that are spread into the atmosphere by the winds, leading to the atmospheric distribution of these elements [3,39,[41][42][43]. Air pollution decreased over the years, even though the drop in the content of elements in Factor 2 was not as clear as in Factor 1 ( Figure 5). In 2002, the enrichment with these elements can be seen along the Vardar River valley, which extends from Skopje to Veles and continues southeastward toward Greece. Due to the frequent flow of air masses in both directions in summer and winter, there was an accelerated diversification of emissions from industry. In 2002, the pollution was mostly connected with the lead and zinc smelter plant in Veles-which had operated until 2003 and with the steel-work smelter in Skopje [2]. Factor 3 (Ni, Cr) contains elements that are usually associated with air pollution, but also influenced by a natural factor, such as the lithological background ( Figure 8). The enrichment of nickel is geogenic, reaching the locality of Groot near the town of Vales in the central part of the country. This location also has Paleogene flysch sediments and Neogene sediments, whose particles can be spread by the wind along the valleys of the Vardar River and Crna River [53]. Similar findings have also been made in previous studies in Macedonia [25,30,81]. The main anthropogenic In 2005, the ER of the elements from Factor 2 is still connected with pollution by the lead and zinc smelter plant in Veles despite its closure in 2003. The high values in 2005 for these elements are also partly associated with the steel-work smelter in Skopje [32]. Leaded gasoline, which was widely used during that period, could also be a possible source of lead in the surface-atmosphere of Macedonia [32]. In contrast with the moss survey in 2002, enrichment with these elements in 2005 was also observed in the area of Bitola in the eastern part of the country. According to the previous study, the Neogene vulcanites were characterized by the natural enrichment of Cd in the region of Berovo-Vladimirovo [42].
In 2010, due to the reactivation of Pb-Zn mines (Sasa, Toranica, Zletovo), the enrichment of Cd and Zn were observed in the eastern part of Macedonia. In the period of 2001/2002 until 2006/2007, the mines were not active, but they were reactivated in 2006 (Sasa) and 2007 (Toranica, Zletovo). In the area of the mines, there were millions of tons of waste material, which can be dispersed into the biosphere by the wind [41][42][43]82]. Furthermore, the enrichment in 2010 was also observed in the area of the steel-work smelter in Skopje. In the surveys from 2010, the contents of Cd and Pb decreased because of the reduced use of leaded gasoline for cars [8].
Factor 3 (Ni, Cr) contains elements that are usually associated with air pollution, but also influenced by a natural factor, such as the lithological background ( Figure 8). The enrichment of nickel is geogenic, reaching the locality of Groot near the town of Vales in the central part of the country. This location also has Paleogene flysch sediments and Neogene sediments, whose particles can be spread by the wind along the valleys of the Vardar River and Crna River [53]. Similar findings have also been made in previous studies in Macedonia [25,30,81]. The main anthropogenic source of these elements is related to the Kavadarci region and includes the ferronickel smelter, the slag dump, and the open pit nickel mine [25,30,31,82]. Generally, the high chromium content belongs to the region near Tetovo, an area of a previous ferrochromium smelter and a slag deposit near the high-melting ferrochromium plant [2].
There were fluctuations in the contents of Ni and Cr, according to the sampling campaign ( Figure 5). The lowest contents were observed in 2002; the highest was in 2005. In samples from moss surveys in 2002, the element contents from the association Factor 3 were related to the geogenic distribution associated with Paleogene flysch sediments and Neogene sediments. The high content was also be observed in the area of Tetovo and Kavadarci, which were associated with a previous ferrochrome smelter and the abandoned ferronickel smelter, respectively. Values for Ni in 2005 were much higher than the values for 2002, which was related to reactivation and increasing production capacity of the ferronickel smelter near Kavadarci in 2004 [30,33,80]. Since 2005, the factory has increased its ore processing capacity, as it also extracts material from Albania, Turkey, and Indonesia [34]. However, the contents of Ni and Cr in the mosses collected in 2010 slightly decreased, which can be linked to a renewed decline in the production capacity at the smelter.
Factor 4 (Ba, K, Sr) is the association of elements that are naturally distributed and whose distribution is not related to industrial and urban activities (Figure 9). The contents of these elements are practically unaltered. The elements are enriched in areas of felsic volcanic rocks (andesite) and their pyroclastic rocks. The contents of these elements slightly increased according to the sampling campaign. In 2002 the enrichment was observed in Kratovo, Ovče Pole, and Kavadarci, where Neogene (Ng) magmatic rock and Paleogene cluster dominates. Another maximum was linked to the southern Macedonia (Kožuf and Mariovo regions) and Strumica area where rocks of Quaternary (Q) sediments dominated. In 2005, the highest values of these elements were found in the area of Kratovo (rocks of Neogene volcanism) [44]; in 2010, the spatial distribution was similar to the distribution from 2002.   Factor 5 (Br, I) includes elements primarily affected by marine influence [74]. Along the Vardar River valley during the winter period, the Jugo wind blows from the northern parts of Africa across the Mediterranean and the Aegean Sea. Consequently, an exponential decrease of Br and I with the distance from the coastline was observed. This phenomenon has already been studied in Norway, where high contents of halogens were observed in soils and mosses corresponding to an exponential decline with distance from the coastline [74,83]. The limestone and dolomite (Paleozoic and Mesozoic carbonates) present in the western part of Macedonia are enriched in iodine, which also explains the elevated values of these elements in the mosses [32]. In 2010, the contents of Br and I were the highest, which means that the wind was stronger in 2010 ( Figure 10).

Conclusions
The moss analysis is a valuable method for monitoring atmospheric deposition of trace elements in Macedonia mainly because it provides a cheap, effective alternative to deposition analysis for identification of areas at risk from high atmospheric deposition fluxes of heavy metals and can play an important role in identifying spatial and temporal trends in atmospheric heavy metal pollution across the country. In this contribution, the results of three air pollution moss surveys of heavy metals are presented and compared. The first systematic study of heavy metals atmospheric pollution using the moss technique was concluded in 2002 [1,2]. This study was followed by two more samplings in 2005 [32] and 2010 [6,22], collecting moss samples from the

Conclusions
The moss analysis is a valuable method for monitoring atmospheric deposition of trace elements in Macedonia mainly because it provides a cheap, effective alternative to deposition analysis for identification of areas at risk from high atmospheric deposition fluxes of heavy metals and can play an important role in identifying spatial and temporal trends in atmospheric heavy metal pollution across the country. In this contribution, the results of three air pollution moss surveys of heavy metals are presented and compared. The first systematic study of heavy metals atmospheric pollution using the moss technique was concluded in 2002 [1,2]. This study was followed by two more samplings in 2005 [32] and 2010 [6,22], collecting moss samples from the same locations. In this study, the results of chemical analyses obtained from the aforementioned studies were statistically processed using bivariate and multivariate statistical techniques. By the comparison of all three moss surveys, it can be concluded that some of the potentially toxic elements ( The largest anthropogenic impact of air pollution with heavy metals was found in the vicinity of the ferronickel smelter near Kavadarci (Ni and Cr), the abandoned lead-zinc smelter near Vales, and the lead and zinc mines in the vicinity of Macedonska Kamenica, Probištip, and Kriva Palanka (Cd, Pb, and Zn). According to factor analysis, three geogenic associations following the lithology and the composition of the soil (Factor 1, which includes Al, As, Co, Cs, Fe, Hf, Na, Rb, Sc, Ta, Th, Ti, U, V, Zr, and RE; Factor 4 with Ba, K, and Sr; Factor 5 with I and Br), one geogenic-anthropogenic association (Factor 3 with Cr and Ni), and one anthropogenic (Factor 2, which includes Cd, Pb, Sb, and Zn) association were obtained.
According to the sampling campaign (2002/2005/2010), the content of elements associated with Factor 1 and Factor 5 changed the most. Factor 1 decreased considerably according to the sampling campaign, while the opposite trend was observed for Factor 5. These two factors show geogenic distribution, mostly related to natural background and weather conditions. The contents of the elements combined in The largest anthropogenic impact of air pollution with potentially toxic metals was the ferronickel smelter plant in the vicinity of Kavadarci (Ni and Cr), as well as three reactivated lead and zinc mines near the towns of Probištip, Makedonska Kamenica and Kriva Palanka (Cd, Pb, and Zn). The distribution of analyzed elements in the Republic of Macedonia was also related to the lithology and the composition of the soil.
This work is essential for monitoring future trends at a high spatial resolution and provides a useful tool for possible modeling the atmospheric deposition fluxes. Environmental monitoring programs, such as the moss survey are also appropriate to regulatory bodies of the country, with the aim to prevent the quality of the environment from deteriorating or ensure that its quality is improved.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1. Matrix of correlation coefficients (n = 216), a group of 16 selected principally natural distributed elements).

Al
As  Values in the range 0.5-0.7 (good association) are underlined, and in the range 0.7-1.0 (strong association) are bolded; Box-Cox transformed values used.