Assessing temporal trends of trace metal concentrations in mosses over France between 1996 and 2011: A flexible and robust method to account for heterogeneous sampling strategies

Environmental Pollution 220 (2017) 828e836

Contents lists available at ScienceDirect

Environmental Pollution journal homepage: www.elsevier.com/locate/envpol

Assessing temporal trends of trace metal concentrations in mosses over France between 1996 and 2011: A flexible and robust method to account for heterogeneous sampling strategies*
a, Aude Pascaud c, d, Emeline Lequy a, *, Nicolas Dubos b, Isabelle Witte c ,d a
bastien Leblond
phane Sauvage , Se Ste a

Natural Heritage Department, National Museum of Natural History, 12 rue Buffon, F-75005, Paris, France Centre d'Ecologie et des Sciences de la Conservation (CESCO UMR 7204) & M
ecanismes adaptatifs et
evolution (MECADEV UMR 7179), Sorbonne Universit
es, MNHN, CNRS, UPMC, CP51, 55 rue Buffon, 75005 Paris, France c Mines Douai, D
epartement Sciences de l'Atmosph ere et G
enie de l’Environment, SAGE, F-59508, Douai, France d Universit
e de Lille, F-59650, Villeneuve-d’Ascq, France b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 16 August 2016 Received in revised form 17 October 2016 Accepted 22 October 2016 Available online 9 November 2016

Air quality biomonitoring has been successfully assessed using mosses for decades in Europe, particularly regarding heavy metals (HM). Assessing robust temporal variations of HM concentrations in mosses requires to better understand to what extent they are affected by the sampling protocol and the moss species. This study used the concentrations of 14 elements measured during four surveys over 15 years in France. Analyses of variance (ANOVA) and a modeling approach were used to decipher temporal variations for each element and adjust them with parameters known to affect concentrations. ANOVA followed by post hoc analyses did not allow to estimate clear trends. A generalized additive mixed modeling approach including the sampling period, the collector and the moss species, plus quadratic effects, was used to analyze temporal variations on repeated sampling sites. This approach highlighted the importance of accounting for non-linear temporal variations in HM, and adjusting for confounding factors such as moss species, species-specific differences between sampling periods, collector and methodological differences in sampling campaigns. For instance, lead concentrations in mosses decreased between 1996 and 2011 following quadratic functions, with faster declines for the most contaminated sites in 1996. On the other hand, other HM showed double trends with U-shaped or hill-shaped curves. The effect of the moss was complex to handle and our results advocate for using one moss species by repeated site to better analyze temporal variations. © 2016 Elsevier Ltd. All rights reserved.

Keywords: Moss survey Bio-monitoring Generalized additive mixed modeling Heavy metals

1. Introduction Air quality has been a major issue tackled by the European authorities for decades, and recently by the Directive 2008/50/EC (EU, 2008), or the Clean Air Policy Package adopted in 2013. Air pollutants are suspected to induce respiratory pathologies especially in the case of particulate matter (Brook et al., 2010; III and Dockery, 2006; WHO, 2014), and they also threaten natural ecosystems and agricultural systems by accumulation of toxic components, notably heavy metals (HM) (Kidd et al., 2015; Nagajyoti et al., 2010;

*

This paper has been recommended for acceptance by Klaus Kummerer. * Corresponding author. E-mail address: [email protected] (E. Lequy).

http://dx.doi.org/10.1016/j.envpol.2016.10.065 0269-7491/© 2016 Elsevier Ltd. All rights reserved.

Sandilyan and Kathiresan, 2014). Biomonitoring studies have been used all over Europe to monitor air quality by measuring HM concentrations in several organisms. The ICP-Vegetation (International Coordinated Program on Effects of Air Pollution on Natural Vegetation and Crops) was created in the late 1980s under the Convention on Long-Range Transboundary Air Pollution. The ICPVegetation uses mosses to bio-monitor air quality in HM, Indeed, mosses have numerous physiological and morphological properties making them efficient bio-monitors of air pollution (Bates, 1992; Harmens et al., 2012; Markert et al., 2003; Ruhling and Tyler, €der et al., 2010). They lack proper root and vascular 1968; Schro systems, and they have no developed cuticle on their leaves. This makes mosses very prone to take up elements from wet and dry deposition, and in addition mosses have a well-developed cation

E. Lequy et al. / Environmental Pollution 220 (2017) 828e836

€rv, 1955). Concentrations in mosses exchange properties (Puustja are higher than in atmospheric deposition and usually above detection levels of analyzers. Besides, their ubiquity and easysampling make it possible to collect mosses on very large territories. Moss bio-monitoring of air quality has already been successfully conducted (Ares et al., 2011; Balabanova et al., 2010; Bekteshi et al., 2015; Berg and Steinnes, 1997; Genoni et al., 2000; Herpin et al., 2001; Varela et al., 2014). For example, a global decrease of HM concentrations in mosses was shown over Europe, suggesting air quality improves (Harmens et al., 2015). Furthermore, relationships between modelled atmospheric deposition and concentrations in mosses were found for lead and cadmium, (Harmens et al., 2012; €der et al., 2014, 2010). Schro However, unlike physical devices used to directly measure atmospheric deposition, mosses are living organisms affected by their environment. Several aspects of the sampling protocol used by the ICP-Vegetation were proven to affect element concentrations in mosses, such as the sampling period, the moss species, or the
ndez et al., 2015; Ferna
ndez and intra-site variability (Ferna Carballeira, 2002; Lequy et al., 2016). In this respect, Dołe˛ gowska and Migaszewski (2015) and Lequy et al. (2016) stressed the need to account for protocol and environmental parameters to better interpret results of moss bio-monitoring. Yet previous studies describing temporal trends over successive campaigns (Harmens et al., 2015) did not include features of the sampling protocol, despite their contribution to the variability of the concentrations in mosses. The objectives of this article are to analyze trends of concentrations in mosses in the French contribution to ICP-Vegetation by (i) using the method usually applied by the ICP-Vegetation members on the data to assess trends and (ii) propose a model that analyze temporal trends of concentrations in mosses and that takes into account the possible effects of the protocol, here the collector, the sampling period and the moss species. To do so, we will (i) use ANOVA and post-hoc Tukey's Honest Significant Differences to globally assess the evolution of concentrations in mosses, and (ii) use mixed modeling on repeated measurements of concentrations during the four French surveys as response variable, with explanatory variables based on features of the protocol. 2. Materials and methods

829

measured since 2000, and As and Sb which were not available in the 2000 survey. In the BRAMM dataset, 280 sites are considered a repetition of measurements over the four surveys (Fig. 1). Such sites were less than 500 m apart between surveys, in similar biotopes (here, the type of forest). Indeed, the BRAMM program only collect mosses under forests, which corresponds to the optimal biotope of the 5 moss species recommended by the ICP-Vegetation guidelines. These 280 sites constitute the “repeated sites” dataset. 2.1.2. Quality of the chemical analyzes The accuracy of the concentration was assessed for the four campaigns using a reference material, following the standardized procedure described in the AFNOR ISO 5472-25 (Lequy et al., 2016). The reference material is based on a large sample of the moss Pleurozium schreberi, named M2, collected in 1993 in a forested background site in Finland as described by Steinnes et al. (1997). For each element, the reference values were computed from the analysis of M2 by 16 European laboratories in 2011. Then the bias between the chemical analyses form a laboratory is determined relatively to the reference value. When the bias is significant, concentrations of the element i were corrected proportionally to the reference concentrations as follows. Corrected concentration_survey_i ¼ measured concentration_survey_i * M2_reference_2011/M2measured by the French laboratory_survey_i. A value by default were attributed to values under the detection or the quantification limits. This value by default was calculated as the minimum of the distribution of the considered BRAMM survey divided by three. Elements As and Sb had no available measurement in 2000. 2.2. Analytical approach and statistics 2.2.1. ANOVA and Tukey's Honest significance difference To highlight an effect of the time on the concentrations in mosses, ANOVA was performed on the four surveys, including all the sampling sites. Differences between the surveys and the regions were assessed using the post hoc test Tukey's Honest Significant Differences (HSD). As all the studied elements were log-normally distributed, they were log-transformed to validate the normality hypothesis needed to use ANOVA and Tukey's HSD. This analysis

2.1. The BRAMM dataset 2.1.1. The BRAMM program The BRAMM program (Bio-monitoring of HM atmospheric deposition using mosses) is the French contribution to the ICPVegetation. It consists in four surveys conducted in 1996, 2000, 2006 and 2011, described in Lequy et al. (2016). They include between 449 and 559 sites. The ICP-Vegetation guidelines recommend to sample 5 moss species as gametophytes: Hylocomium splendens (Hedw.) Schimp., Hypnum cupressiforme Hedw., Pseudoscleropodium purum (Hedw.) M.Fleisch. in Engl. et al., Pleurozium schreberi (Brid.) Mitt., and Thuidium tamariscinum (Hedw.) Schimp., respectively referred to as Hs, Hc, Pp, Ps, and Tt. Details about the sampling protocols and analytical procedures are described in (Harmens, 2010). Different other data were recorded, notably the site identification, the year of the survey, the date of sampling, the moss species, and the name of the collector of each sampling site for each survey. This study focuses on the following HM and other elements: Al, As, Ca, Cd, Cr, Cu, Fe, Hg, Na, Ni, Pb, Sb, V and Zn. These elements were measured in the 4 surveys, except Hg which has been

Fig. 1. Distribution of the repeated sites of the BRAMM datasets relatively to all the sampling sites pooled from the four surveys.

830

E. Lequy et al. / Environmental Pollution 220 (2017) 828e836

Table 1 Summary of the variables used by the mixed modeling approach. Df stands for degrees of freedom, and doy for “day of the year”. Name

Variables

Expression of the fixed effects in the model

Df

Simplest model (M-t)

Fixed effect: Time

Time

4

Time*class

13

(Time þ Time2)*class

18

Model including the classes, quadratic effects, Fixed effects: Same as M-tcq, plus moss species, interaction moss species and and the protocol time, and day of the year (spline) (M-tcqp) Random effects: same as M-t, plus collector

(Time þ Time2)*class þ Time*moss species þ s(doy)

29

Model including the classes, the protocol, and Fixed effects: Same as R-tcqp with an interaction between sampling period (spline) and moss species. quadratic effects Random effects: same as M-tcqp (M-tcqpþ)

(Time þ Time2)*class þ Time*moss species þ s(doy*moss species)

37

Random effect: site, year of survey Model including the classes (M-tc)

Fixed effects: Same as M-t, plus classes, and interaction classes and Time Random effect: same as M-t

Model including the classes and quadratic effects (M-tcq)

Fixed effect: Same as M-tc, plus quadratic effect on the interaction classes and Time Random effect: same as M-t

allows to establish differences between surveys, but cannot include other covariates such as the protocol features described in the former section. Normality of the residuals and homoscedasticity were validated after visual inspection. These analyses were performed using the R freeware version 3.21.1 (R Core Team, 2015). 2.2.2. Mixed modeling Repeated measurement every 4e5 years on sampling sites induces a structure in the data that violates the assumption of independence between samples. For this reason we used mixed models, which allow for dependence between samples (Zuur et al., 2009). Data was log-transformed to meet the demand of normality and homoscedasticity of the residuals. A series of models was applied to (i) select the model best describing the temporal trends and (ii) assess the effect of the features of the protocol on the models’ results (Table 1). We used Generalized Additive Mixed Models (GAMM) from the MGCV package (Wood, 2006) implemented in R version 3.2.1 (R Core Team, 2015). This enabled us to describe the relationship between the response variable (i.e., the element concentration) and Time, while accounting for non-linear confounding factors. We treated each element separately. However, all moss species were analyzed simultaneously, so that the modeling of adjustment variables relied on the maximum amount of data. 2.2.2.1. Fixed effects. Time: Time (in years) was represented by a numeric vector of 0, 4, 10, and 15 years after the 1996 survey and included as fixed effect. Classes of concentration: for each sample, the sampling period was the day of year during which mosses were collected. From the “repeated sites” dataset, classes of 20th percentiles (five quartiles, with lowest and largest concentrations in the 1st and 5th classes, respectively) were computed on the basis of element concentrations in mosses in the 1996 survey (or in the 2000 survey for Hg). These classes were attributed to all the repeated sampling sites over metropolitan France. These classes very likely represent a gradient of industrial activity (As, Cd, Cr, Cu, Hg, Ni, Pb, Sb, and Zn) or deposits of dust outburst from Africa or local erosion (Al, Fe, Ca), or a gradient of sea spray (Na), leading to different patterns of element concentrations. These classes present different climatologies and various source types. Interaction between classes of concentration and Time: in order

to better understand the temporal trends relatively to their initial concentration, we included an interaction term between Time and concentration class. This enabled to estimate one coefficient for each concentration class individually. Quadratic effect: in case the temporal variation in element concentration was non-linear, for each element, we built a second model that included a quadratic effect of Time. 2.2.2.2. Fixed effect related to the sampling protocol. We included a number of adjustment variables to our models, to control for confounding effects related to the sampling protocol. Day of the year (doy): element concentrations may vary nonlinearly within a year, and these variations may differ between species. We therefore included a “day of the year” effect, with a spline function fitted for each species. The “day of year” variable needed to be scaled for the models to converge. Moss species: we controlled for between-species differences in element concentrations by including a moss species fixed term. Interaction between moss species and time: because temporal trends in element concentrations may vary between species, we included an interaction term between species and Time. 2.2.2.3. Random variables. To account for changes in sampling procedures between surveys that were not related to the temporal trend, we included “year of survey” as random variable. We finally accounted for between sites, and between moss collector differences by adding a site and a collector random effect, respectively. Before using mixed models, the possible co-variables (moss species, sampling period) were checked to have significant effects on element concentrations by Kruskal-Wallis. The tests were significant for the 14 elements considered in this study, hence sampling period and moss species were used as co-variables in the models. Models were compared by ANOVA using the MGCV package, based on AIC and generalized likelihood test, to select the best model for each element. Models were compared on the basis of Akaike Information Criterion (AIC), the lowest AIC being attributed to models with the highest statistical support (Burnham and Anderson, 2004). Visual inspection of residual plots did not reveal any major deviations from homoscedasticity or normality. Effects were considered significant for p-values <5%.

E. Lequy et al. / Environmental Pollution 220 (2017) 828e836

3. Results

831

Table 3 ANOVA between the 5 models for the 14 elements. AIC is the Akaike Index Criterion. P-values in bold are <5%. AIC values in bold are the lowest for each element and allow to select the best fitting model for each element. N ¼ 280 sites.

3.1. Temporal variations of the 14 elements analyzed in the French surveys by ANOVA and Tukey One factor ANOVA of the log-transformed element concentrations with the year were significant for all the elements analyzed in this study. This highlights strong effects of the surveys. The Tukey's HSD test highlighted differences of element concentrations between the surveys and made it possible to establish a decreasing trend for Na and Pb over France (Table 2), but results were more intricate for the other elements. The elements Sb and V seemed to decrease, while Cu and Ni seemed to increase, and the elements Al, As, Ca, Cd, Cr, Hg, and Zn did not exhibit any clear trend.

Element

Model

AIC

p-value

Element

Model

AIC

p-value

Al

M-t M-tc M-tcq M-tcqp M-tcqpþ

3763 2147 1974 1912 1911

NA 0.0Eþ00 0.0Eþ00 2.9E-13 2.5E-02

Hg

M-t M-tc M-tcq M-tcqp M-tcqpþ

2596 120 21 75 64

NA 0.0Eþ00 0.0Eþ00 8.6E-12 8.4E-01

As

M-t M-tc M-tcq M-tcqp M-tcqpþ

2511 2063 1917 1920 1934

NA 0.0Eþ00 0.0Eþ00 6.5E-02 9.6E-01

Na

M-t M-tc M-tcq M-tcqp M-tcqpþ

3245 1825 1816 1559 1551

NA 0.0Eþ00 1.4E-03 0.0Eþ00 2.5E-03

Ca

M-t M-tc M-tcq M-tcqp M-tcqpþ

2725 561 542 388 398

NA 0.0Eþ00 2.4E-05 0.0Eþ00 6.7E-01

Ni

M-t M-tc M-tcq M-tcqp M-tcqpþ

3117 2305 1651 1559 1558

NA 0.0Eþ00 0.0Eþ00 0.0Eþ00 2.7E-02

Cd

M-t M-tc M-tcq M-tcqp M-tcqpþ

2719 1643 1522 1420 1440

NA 0.0Eþ00 0.0Eþ00 0.0Eþ00 8.6E-01

Pb

M-t M-tc M-tcq M-tcqp M-tcqpþ

2273 1096 1019 869 881

NA 0.0Eþ00 0.0Eþ00 0.0Eþ00 8.3E-01

Cr

M-t M-tc M-tcq M-tcqp M-tcqpþ

2314 1796 1703 1658 1664

NA 0.0Eþ00 0.0Eþ00 3.4E-10 3.2E-01

Sb

M-t M-tc M-tcq M-tcqp M-tcqpþ

2547 1556 994 977 987

NA 0.0Eþ00 0.0Eþ00 3.8E-05 6.7E-01

Cu

M-t M-tc M-tcq M-tcqp M-tcqpþ

1835 455 408 367 361

NA 0.0Eþ00 4.7E-11 2.3E-09 6.5E-03

V

M-t M-tc M-tcq M-tcqp M-tcqpþ

2272 1520 1296 1263 1261

NA 0.0Eþ00 0.0Eþ00 7.1E-08 2.1E-02

Fe

M-t M-tc M-tcq M-tcqp M-tcqpþ

3499 1864 1780 1754 1757

NA 0.0Eþ00 0.0Eþ00 1.4E-06 1.1E-01

Zn

M-t M-tc M-tcq M-tcqp M-tcqpþ

2584 1043 1005 819 848

NA 0.0Eþ00 4.6E-09 0.0Eþ00 1.0E-01

3.2. Mixed models and effect of the sampling protocol and the concentration classes on temporal variations This section uses the 280 ‘repeated sites’ dataset. 3.2.1. Effects of the initial concentrations and of the features of the protocol on the quality of the models: selection of the best fitting model for each element The five models presented were compared from the simplest to the most complex by ANOVA for the 14 elements studied in this article (Table 3). For all elements, the inclusion of classes (M-tc) and quadratic effects (M-tcq) systematically improved model fits. AIC also decreased after adding the protocol features (M-tcqp) for all the elements except for As. The interaction between the sampling period (M-tcqpþ) and the moss species significantly improved model fit for Al, Cu, Na, Ni, and V. Henceforth, for each studied element, the data are better explained by adding the classes of initial concentrations (1996 survey), and adding features of the protocol also improved the quality of the model and the explanation of data variability. 3.2.2. Results of the mixed modeling approach and description of the trends for the 14 elements Following the previous section, the temporal trends were described according to the models including only Time and classes for As, including the protocol without any interaction between the sampling period and the moss species (Mtcqp) for Ca, Cd, Cr, Fe, Hg, Pb, Sb, and Zn, or with this interaction (Mtcqpþ) for Al, Cu, Na, Ni and V (Table 4). The elements showed significantly different concentrations according to the “class” variable. The interaction variables “classes:Time” were significant for the 5 classes for As, Cd, Na, Pb, and Sb, and for between 3 and 4 classes for Al, Ca, Cr, Cu, Fe, Hg, Ni, V, and Zn. For Na and Pb, the “class” interaction with Time indicated that sites with the highest concentrations in element showed the fastest decrease in element. The elements Al, As, Cr, Fe, Hg, and V presented significantly negative and positive coefficients for the variables ‘Interaction classes:Time’ and “Quadratic Interaction

classes:Time”, respectively, for at least the 3 upper classes. This suggests hump-shaped variations for such elements, with an inflexion after 2000 (Fig. 2). On the contrary, Sb presented significantly positive and negative coefficients for the variables ‘Interaction classes:Time’ and “Quadratic Interaction classes:Time”, respectively, for at least 3 classes. This suggests U-shaped variations for Sb (Fig. 2). The Ca concentrations significantly increased in mosses throughout the four surveys for the 3 upper classes of sites. The elements Cd, Cu, Ni, and Zn showed different behaviors of the classes of site 1 and 5. The coefficients of Table 4 and Fig. 2 suggest that the class 1 follows a quite linearly increasing curve, while the class 5 follows a decline until the middle of the 2000s, then an increase until 2011 (see Fig. 2). At least one moss species showed a significant base concentration for Al, Ca, Cd, Cr, Cu, Fe, Pb, V, and Zn, relatively to the Hc moss species, indicating different base levels of each concerned element

Table 2 Survey means of concentrations in mosses for 14 elements (mg.g1) from 1996 to 2011 analyzed by Tukey's HSD. NA stands for not available, i.e. elements for which concentrations could not be measured. Letters a, b, and c show the significantly different values for each element. N ¼ 2025. Survey

Al

As

Ca

Cd

Cr

Cu

Fe

Hg

Na

Ni

Pb

Sb

V

Zn

1996 2000 2006 2011

646b 679b 966a 484c

0.3b NA 0.4a 0.2c

3896b 4657a 4148b 4985a

0.22a 0.17b 0.12c 0.18b

1.7b 1.6b 2.4a 1.5b

5.5b 5.5b 5.7b 6.4a

588b 612b 749a 503c

NA 0.07b 0.09a 0.07b

215a 150b 138b 92c

2.0b 1.9b 2.2a 2.3a

9.8a 5.6b 4.3c 3.7d

0.22a NA 0.21a 0.15b

2.6a 2.6a 2.6a 1.3b

34b 44a 31c 34b

832

E. Lequy et al. / Environmental Pollution 220 (2017) 828e836

Table 4 Coefficient estimates of Time, interaction of time and classes of concentrations (from class 1 to 5), moss species, and interaction of time and moss species, computed by best fitting model for each element. Values in bold were significant (p-value <5%). The values of the concentrations in the variables ‘class’ and ‘moss’ are in log due to the logtransformation of the response variable. The ‘moss’ and ‘moss:Time’ variables were computed relatively to Hc and mossHc:Time, respectively. N ¼ 280 sampling sites. Explanatory Variable

Al

As

Ca

Cd

Cr

Cu

Fe

Hg

Na

Ni

Pb

Sb

V

Zn

Model

Mtcqpþ

Mtcq

Mtcqp

Mtcqp

Mtcqp

Mtcqpþ

Mtcqp

Mtcqp

Mtcqpþ

Mtcqpþ

Mtcqp

Mtcqp

Mtcqpþ

Mtcqp

class class class class class

1 2 3 4 5

5.56 5.92 6.28 6.63 7.31

¡1.93 ¡1.46 ¡1.23 ¡0.89 0.09

7.59 8.04 8.24 8.50 8.93

¡2.72 ¡1.87 ¡1.63 ¡1.43 ¡0.88

¡0.22 0.13 0.45 0.74 1.36

1.05 1.42 1.64 1.86 2.22

5.62 5.95 6.21 6.63 7.20

¡4.12 ¡3.52 ¡3.33 ¡2.99 ¡2.32

4.74 5.11 5.33 5.63 6.33

¡2.96 ¡0.29 0.38 0.78 1.42

1.66 2.00 2.24 2.48 3.00

¡2.36 ¡1.91 ¡1.66 ¡1.35 ¡0.73

0.22 0.59 0.81 1.11 1.64

2.93 3.28 3.44 3.62 4.03

class class class class class

1:Time 2:Time 3:Time 4:Time 5:Time

0.21 0.11 0.12 0.06 0.03

0.24 0.15 0.12 0.16 0.10

0.13 0.07 0.06 0.02 0.01

0.09 ¡0.05 ¡0.08 ¡0.09 ¡0.15

0.19 0.14 0.07 0.04 0.01

0.10 0.03 0.02 0.02 ¡0.07

0.14 0.07 0.09 0.02 0.03

0.31 0.20 0.19 0.14 0.06

¡0.06 ¡0.11 ¡0.15 ¡0.16 ¡0.25

0.86 0.21 0.06 0.04 ¡0.06

0.02 ¡0.06 ¡0.10 ¡0.10 ¡0.17

¡0.19 ¡0.27 ¡0.28 ¡0.26 ¡0.33

0.17 0.10 0.06 0.04 0.01

0.11 0.05 0.03 0.00 ¡0.04

class class class class class

1:Time2 2:Time2 3:Time2 4:Time2 5:Time2

¡0.01 ¡0.01 ¡0.01 ¡0.01 ¡0.01

¡0.02 ¡0.01 ¡0.01 ¡0.01 ¡0.01

0.00 0.00 0.00 0.00 0.00

0.00 0.00 0.01 0.01 0.01

¡0.01 ¡0.01 0.00 0.00 0.00

0.00 0.00 0.00 0.00 0.00

¡0.01 0.00 ¡0.01 0.00 0.00

¡0.01 ¡0.01 ¡0.01 ¡0.01 0.00

0.00 0.00 0.00 0.00 0.01

¡0.04 ¡0.01 0.00 0.00 0.00

0.00 0.00 0.00 0.00 0.01

0.01 0.02 0.02 0.02 0.02

¡0.01 ¡0.01 ¡0.01 ¡0.01 0.00

¡0.01 0.00 0.00 0.00 0.00

mossHs mossPp mossPs mossTt

0.29 0.36 0.11 0.12

NA NA NA NA

¡0.32 0.04 ¡0.09 0.03

0.15 0.05 0.12 0.15

0.09 0.22 0.15 0.14

0.10 0.04 0.10 0.06

0.01 0.19 0.15 0.15

0.03 0.04 0.08 0.00

0.01 0.20 0.17 0.05

0.15 0.01 0.08 0.07

0.05 0.05 0.01 ¡0.13

0.07 0.02 0.05 0.00

0.16 0.30 0.12 0.06

0.03 0.02 0.01 0.29

mossHs:Time mossPp:Time mossPs:Time mossTt:Time

0.01 ¡0.05 0.01 0.02

NA NA NA NA

0.03 0.01 0.01 ¡0.02

¡0.04 ¡0.02 0.02 ¡0.04

0.00 ¡0.02 0.01 ¡0.02

0.00 0.00 0.02 0.02

0.02 ¡0.03 0.01 0.02

0.00 ¡0.02 ¡0.02 0.00

0.03 0.03 0.04 0.06

0.03 ¡0.02 0.00 0.01

¡0.02 ¡0.03 ¡0.04 ¡0.02

0.00 ¡0.02 0.00 ¡0.02

0.00 ¡0.05 0.02 0.02

0.01 0.02 0.01 ¡0.02

in moss species. The interaction moss:Time was significant for all elements except As at least for one moss species relatively to the Hc moss species. Finally, a non-linear effect “period of the year” was found significant for all the elements except As and Sb (Table 5). The effect of the “period of the year” was different between moss species for Cu, Na, Ni, and V. For Ca, Cd, Cr, Fe, Hg, Pb, and Zn, smooth terms showed between 4 and 7 degrees of freedom, indicating between 4 and 7 changes of the curve inflection throughout a year, as illustrated by Fig. 3 for Al, Ca, Na, and Zn.

4. Discussion 4.1. Effect of protocol parameters on the description of temporal trends 4.1.1. Effect of the protocol parameters on the quality of the models Estimations of temporal trends in element concentrations differed between the ANOVA (Table 2) and the generalized additive mixed modeling approach in a majority of cases. The results of the modeling approach indicated that, for each element, temporal trends are function of the classes of sites (depending on the first survey). For this reason, we think that temporal trends assessed without accounting for such classes often mask a more complex reality. The model selection clearly showed an improvement of the model with the inclusion of the protocol information, henceforth significantly improving the description of the element concentrations through time. Previous articles dealing with temporal changes of concentrations in mosses did not include any parameter of the protocol (Barandovski et al., 2015, 2012; Cao et al., 2008;  Harmens et al., 2015, 2010; Nickel et al., 2014; Sakalys et al., 2009). Our study showed the importance of accounting for the protocol to disentangle significant temporal trends and provide robust results.

For this reason, we strongly advise to adjust concentrations with every available and relevant parameter of the protocol when analyzing trends. Besides, the moss species is a very important parameter to take into account and it is of primary importance to deal with it when analyzing trends.

4.1.2. Influence of the quadratic effects on the assessment of temporal trends Annual averages in France between 1996 and 2011 (Fig. 2) are very similar to those found in Europe (Harmens et al., 2015) for Al, Cd, Cu, Fe, Hg, and Pb. A decline was found to occur for Pb at the European scale and in the present study. However, the declines found in the European scale were not as clear for Cd or Hg in the present study. The elements Cd and Hg seem also to show increasing concentrations, at least up to the 2006 survey in France. Besides, Cu concentrations mainly increased in the classes 1 to 4 in France, while a decline was assessed at the European scale. This can be due to a national effect, as France is a country less contaminated by trace metals than Eastern countries (Harmens et al., 2015), so global trends do not necessary correspond to those in France. The main difference with the European temporal trends lay in significant quadratic effects. They indicate non-linear trends for all the elements. This makes it difficult to estimate simple decline or increase between two temporal limits, except for the elements showing clear trends on Fig. 2 (Na and Pb). This would completely mask inflections of the trends, as seen for Al or Sb. This nonlinearity is also acknowledged by the EMEP Program with the use of the Mann-Kendall method to describe non monotonous temporal trends of atmospheric pollutants (Sauvage et al., 2009). This is crucial for scientists and policy-makers, in order to explain inflections of element concentrations in mosses, their causes, and how to possibly handle them to durably reduce HM in the atmosphere and in the ecosystems. The results of this study showed for the first time that for Pb, the most concentrated sites showed the quickest decline. This is a novel

E. Lequy et al. / Environmental Pollution 220 (2017) 828e836

833

Fig. 2. Boxplots of the distributions (0.25, 0.5 and 0.75 quartiles plus outliers after 1.5 interquartile interval) of the concentrations in mosses for 14 elements (mg.g1 DW) after removal of values larger than the median plus three standard deviations for the 14 elements studied in the present article for the four BRAMM surveys (n ¼ 280 sampling sites), except 1996 and 2000 for Hg, and As and Sb, respectively. Points are annual means. The colored curves represent the predicted values obtained from GAMM for the five concentration classes after back-transformation of the modelled data.

834

E. Lequy et al. / Environmental Pollution 220 (2017) 828e836

Table 5 Results of the spline variable “day of the year” (or doy), with or without any interaction with the moss species in function of the best fitting model for each element. df stands for degrees of freedom and represents the number of inflexions of the curve between element concentrations and doy. Bold values were significant (p-value <5%). N ¼ 280 sampling sites. Element

Model

Al

Mtcqp+

As Ca Cd Cr

Mtcq Mtcqp Mtcqp Mtcqp

Cu

Mtcqp+

Fe Hg

Mtcqp Mtcqp

Variable

df

p.value

doy:Hc doy:Hs doy:Pp doy:Ps doy:Tt NA doy doy doy doy:Hc doy:Hs doy:Pp doy:Ps doy:Tt doy doy

5 1 1 1 1 NA 4 4 6 1 1 2 1 4 5 1

5.2E-11 1.2E-01 2.5E-04 1.0E+00 2.6E-03 NA 2.4E-17 3.6E-10 1.2E-10 2.0E-01 8.7E-01 8.4E-03 1.2E-01 1.4E-04 3.8E-06 1.9E-02

and interesting result that illustrates the efficiency of abatement policies applied in France since the middle 1990s. The dynamics of other elements are also very intriguing, as the trends were double for most of elements (increasing then decreasing, or the contrary).

4.2. Recommendations for future analyses of trends The present study confirmed the effects of moss sampling protocol on the assessment of temporal trends for HM concentrations.

Element

Model

Na

Mtcqp+

Ni

Mtcqp+

Pb Sb

Mtcqp Mtcqp

V

Mtcqp+

Zn

Mtcqp

Variable

df

p.value

doy:Hc doy:Hs doy:Pp doy:Ps doy:Tt doy:Hc doy:Hs doy:Pp doy:Ps doy:Tt doy doy doy:Hc doy:Hs doy:Pp doy:Ps doy:Tt doy

8 1 7 1 1 1 1 3 3 1 1 1 4 1 4 1 1 7

1.8E-16 6.9E-01 2.9E-08 9.4E-01 4.2E-02 3.7E-01 3.3E-02 9.9E-05 6.2E-05 3.8E-01 4.6E-02 8.6E-01 5.4E-06 2.1E-01 7.8E-03 9.0E-01 3.7E-02 4.4E-21

Indeed, the present study highlights that the effect of the moss species is very delicate to isolate and they can have contrasted effects on the description of the temporal trends for many of the elements studied in this article Our results also indicated strong effects of the sampling period on the concentrations. For this reason we recommend to sample as much as possible the same moss species during the shortest possible and similar sampling periods on the same sites from a survey to another in order to assess as robust as possible temporal trends on the long term.

Fig. 3. Representation of the spline effect of the sampling period (“day of year” or doy) on the concentrations of elements (here Al and Na with the interaction with Hc moss species), and Ca and Zn (without any interaction with moss species). Sampling period was computed as day of year (doy) that was scaled for the models to converge. The plain and dotted lines represents the average effect of doy on the concentrations in mosses and its standard deviation, respectively.

E. Lequy et al. / Environmental Pollution 220 (2017) 828e836

The BRAMM surveys do not include any growth or biological activity data. As the sampling period had been shown to affect
ndez et al., 2015; Lequy et al., concentrations in mosses (Ferna 2016), we used this variable that can account for other variables, such as, for example, biological activity. Besides, growth data is very delicate to obtain and associated with large uncertainty, as mosses do not grow as vascular plants (Boquete et al., 2014; Leblond et al., 2004). Indeed, as mosses grow on their apex, their bottom decay. In addition, mosses are sampled at one particular date by survey, so punctual growth rate would be difficult to include in the model. Regarding the moss species, those recommended by the ICPVegetation present very different morphologies, such as the size or shape of their leaves, which can affect the contact surface with atmosphere and thus heavy metal bioaccumulation. For this reason, accounting for the moss species is mandatory in the modeling approach. Since mosses were only sampled as gametophytes and since the moss species recommended by the ICP-Vegetation have mainly an asexual reproduction, reproduction cannot affect our results, so it was not accounted for in the models. Real time deposition of heavy metals and other elements is also largely unavailable, especially on 280 sampling sites in France. For this reason, many countries use moss bio-monitoring to make up for the lack of fine grid monitoring of heavy metal deposition. Data about moisture and temperature were not included in the dataset, so they could not be included in the modeling approach. Such variables can affect concentrations in mosses, so, when available, they could be included in the modeling approach. Biological features, such as the leaf morphology, physiology, or growth rate, are still to be studied to refine the modeling approach proposed in this study. In any case, the objective of this study was not as much modeling the biological behavior of moss towards heavy metals as assessing temporal variations of the concentrations measured in bio-monitoring mosses, while accounting for the main potential confounder variables of the sampling protocol. We also strongly advise to analyze trends while accounting for non-linear trends. When possible, we also suggest to use the repeated sites, even though the definition of a repeated site is less strict according to the ICP-Vegetation. We strongly advise to consider a site as a repetition only when it is in the same biotope than that of the previous survey, as now recommended by the ICPVegetation guidelines (Frontasyeva and Harmens, 2014). We suggest to analyze trends not only geographically, but by classes of concentrations, which are usually geographically distributed, and less arbitrary than man-made borders while HM are associated with long-range transboundary pollution. This would validate the method proposed in this study at the European scale and allow to compare the results with previous analyses. 5. Conclusion This study is the first attempt to analyze temporal variations by comparing an ANOVA approach and a generalized additive mixed modeling approach. The models included variables of the protocol (sampling period, collector, and moss species) and allowed for nonlinear temporal variations. The use of ANOVAs in detecting temporal trends of element concentrations in mosses may be over simplistic in a majority of cases. This approach cannot handle non-linear variations, which may lead to spurious results in most studies. In the case of multiple surveys with inconstant sampling strategies, a modeling approach enables to adjust for confounding effects. We recommend the use of mixed models, especially GAMMs, to enable accounting for nonlinear variations in element such as within-year variations, and collecting different moss species by different collectors. In this paper, the 14 analyzed elements showed non-linear

835

temporal variations with a global decline for Pb and at least twotrends following quadratic functions between 1996 and 2011. Classes of concentrations are easy to compute and including them in the models significantly improved the quality of fit for the 14 elements. Besides, all elements, except As, were significantly affected by the variables of the protocol, which also significantly improved the quality of fit of the corresponding element. Therefore, temporal trends were better described. This advocates for accounting for such variables when analyzing temporal variations, and if possible, reduce the sampling period to better isolate longterm temporal variations. The moss species was delicate to isolate. This also advocates for collecting one moss species by repeated sampling site. Temporal trends assessed at large spatial scale could benefit from our results. Our modeling approach was statistically validated with only simple or quickly computed variables. Therefore, our modeling approach is robust, and can be adapted to various areas. We propose this method to assess temporal variations in other countries or at the European scale provided (i) measurements were carried out on repeated sampling sites and (ii) some basic information about the protocol are included in databases. Acknowledgments This study was funded by the ADEME project no1462C0012 (French Agency for Environment and Energy Management). We are very grateful to all the collectors who contributed to the BRAMM database during the four surveys, and to those who took care of the samples in the laboratory. We are also in debt to the former mans Letrouit, Catherine agers of the BRAMM program, Marie-Agne s and Rausch, Sandrine Gombert-Courvoisier, Laurence Galsomie Xavier Laffray. We also want to thank the anonymous reviewers who helped us improve the manuscript. References III, C.A.P., Dockery, D.W., 2006. Health effects of fine particulate air pollution: lines that connect. J. Air Waste Manag. Assoc. 56, 709e742. http://dx.doi.org/10.1080/ 10473289.2006.10464485.  Balabanova, B., Stafilov, T., Ba ceva, K., Sajn, R., 2010. Biomonitoring of atmospheric pollution with heavy metals in the copper mine vicinity located near Radovis, Republic of Macedonia. J. Environ. Sci. Health Part A 45, 1504e1518. http:// dx.doi.org/10.1080/10934529.2010.506097.  Barandovski, L., Frontasyeva, M.V., Stafilov, T., Sajn, R., Pavlov, S., Enimiteva, V., 2012. Trends of atmospheric deposition of trace elements in Macedonia studied by the moss biomonitoring technique. J. Environ. Sci. Health Part A 47, 2000e2015. http://dx.doi.org/10.1080/10934529.2012.695267.  Barandovski, L., Frontasyeva, M.V., Stafilov, T., Sajn, R., Ostrovnaya, T.M., 2015. Multielement atmospheric deposition in Macedonia studied by the moss biomonitoring technique. Environ. Sci. Pollut. Res. 22, 16077e16097. http:// dx.doi.org/10.1007/s11356-015-4787-x. Bates, J.W., 1992. Mineral nutrient acquisition and retention by bryophytes. J. Bryol. 17, 223e240. http://dx.doi.org/10.1179/jbr.1992.17.2.223. Bekteshi, L., Lazo, P., Qarri, F., Stafilov, T., 2015. Application of the normalization process in the survey of atmospheric deposition of heavy metals in Albania through moss biomonitoring. Ecol. Indic. 56, 50e59. http://dx.doi.org/10.1016/ j.ecolind.2015.03.001. Berg, T., Steinnes, E., 1997. Use of mosses (Hylocomium splendens and Pleurozium schreberi) as biomonitors of heavy metal deposition: from relative to absolute deposition values. Environ. Pollut. 98, 61e71. http://dx.doi.org/10.1016/S02697491(97)00103-6.
ndez, J.A., 2014. Effect of age on the Boquete, M.T., Aboal, J.R., Carballeira, A., Ferna heavy metal concentration in segments of Pseudoscleropodium purum and the biomonitoring of atmospheric deposition of metals. Atmos. Environ. 86, 28e34. http://dx.doi.org/10.1016/j.atmosenv.2013.12.039. Brook, R.D., Rajagopalan, S., Pope, C.A., Brook, J.R., Bhatnagar, A., Diez-Roux, A.V., Holguin, F., Hong, Y., Luepker, R.V., Mittleman, M.A., Peters, A., Siscovick, D., Smith, S.C., Whitsel, L., Kaufman, J.D., 2010. Particulate Matter Air Pollution and Cardiovascular Disease An Update to the Scientific Statement From the American Heart Association. Circulation 121, 2331e2378. doi:10.1161/ CIR.0b013e3181dbece1. Burnham, K.P., Anderson, D.R. (Eds.), 2004. Model Selection and Multimodel Inference. Springer New York, New York, NY. Cao, T., An, L., Wang, M., Lou, Y., Yu, Y., Wu, J., Zhu, Z., Qing, Y., Glime, J., 2008. Spatial

836

E. Lequy et al. / Environmental Pollution 220 (2017) 828e836

and temporal changes of heavy metal concentrations in mosses and its indication to the environments in the past 40 years in the city of Shanghai. China. Atmos. Environ. 42, 5390e5402. http://dx.doi.org/10.1016/ j.atmosenv.2008.02.052. Dołe˛ gowska, S., Migaszewski, Z.M., 2015. Plant sampling uncertainty: a critical review based on moss studies. Environ. Rev. 23, 151e160. http://dx.doi.org/ 10.1139/er-2014-0052. EU, 2008. Directive 2008/50/EC of the European Parliament and of the Council of 21 May 2008 on Ambient Air Quality and Cleaner Air for Europe.
ndez, J.A., Carballeira, A., 2002. Biomonitoring metal deposition in Galicia (NW Ferna Spain) with mosses: factors affecting bioconcentration. Chemosphere 46, 535e542. http://dx.doi.org/10.1016/S0045-6535(01)00060-1.
ndez, J.A., Boquete, M.T., Carballeira, A., Aboal, J.R., 2015. A critical review of Ferna protocols for moss biomonitoring of atmospheric deposition: sampling and sample preparation. Sci. Total Environ. 517, 132e150. http://dx.doi.org/10.1016/ j.scitotenv.2015.02.050. Frontasyeva, M., Harmens, H., The participants of the ICP Vegetation, 2014. Heavy Metals, Nitrogen and POPs in European Mosses: 2015 Survey. Genoni, P., Parco, V., Santagostino, A., 2000. Metal biomonitoring with mosses in the surroundings of an oil-fired power plant in Italy. Chemosphere 41, 729e733. http://dx.doi.org/10.1016/S0045-6535(99)00457-9. Harmens, H., 2010. Monitoring of atmospheric deposition of heavy metals, nitrogen and POPs in Europe using Bryophytes (Monitoring manual). Harmens, H., Norris, D.A., Steinnes, E., Kubin, E., Piispanen, J., Alber, R., Aleksiayenak, Y., Blum, O., Cos¸kun, M., Dam, M., De Temmerman, L.,
ndez, J.A., Frolova, M., Frontasyeva, M., Gonz
Ferna alez-Miqueo, L.,
ska, K., Jeran, Z., Korzekwa, S., Krmar, M., Kvietkus, K., Leblond, S., Grodzin  kovsk
Liiv, S., Magnússon, S.H., Man a, B., Pesch, R., Rühling, Å., Santamaria, J.M., € der, W., Spiric, Z., Suchara, I., Tho €ni, L., Urumov, V., Yurukova, L., Schro Zechmeister, H.G., 2010. Mosses as biomonitors of atmospheric heavy metal deposition: spatial patterns and temporal trends in Europe. Environ. Pollut. 158, 3144e3156. http://dx.doi.org/10.1016/j.envpol.2010.06.039. Harmens, H., Ilyin, I., Mills, G., Aboal, J.R., Alber, R., Blum, O., Cos¸kun, M., De
Figueira, R., Frontasyeva, M., Godzik, B.,
ndez, J.A., Temmerman, L., Ferna Goltsova, N., Jeran, Z., Korzekwa, S., Kubin, E., Kvietkus, K., Leblond, S., Liiv, S.,  kovsk
Magnússon, S.H., Man a, B., Nikodemus, O., Pesch, R., Poikolainen, J., €der, W., Spiric, Z., Stafilov, T., Radnovi
c, D., Rühling, Å., Santamaria, J.M., Schro €ni, L., Turcsa
nyi, G., Yurukova, L., Steinnes, E., Suchara, I., Tabors, G., Tho Zechmeister, H.G., 2012. Country-specific correlations across Europe between modelled atmospheric cadmium and lead deposition and concentrations in mosses. Environ. Pollut. 166, 1e9. http://dx.doi.org/10.1016/ j.envpol.2012.02.013. Harmens, H., Norris, D.A., Sharps, K., Mills, G., Alber, R., Aleksiayenak, Y., Blum, O., Cucu-Man, S.-M., Dam, M., De Temmerman, L., Ene, A., Fern
andez, J.A., Martinez-Abaigar, J., Frontasyeva, M., Godzik, B., Jeran, Z., Lazo, P., Leblond, S., Liiv, S.,  kovska
, B., Karlsson, G.P., Piispanen, J., Poikolainen, J., Magnússon, S.H., Man Santamaria, J.M., Skudnik, M., Spiric, Z., Stafilov, T., Steinnes, E., Stihi, C., €ni, L., Todoran, R., Yurukova, L., Zechmeister, H.G., 2015. Heavy Suchara, I., Tho metal and nitrogen concentrations in mosses are declining across Europe whilst some “hotspots” remain in 2010. Environ. Pollut. 200, 93e104. http:// dx.doi.org/10.1016/j.envpol.2015.01.036. Herpin, U., Siewers, U., Kreimes, K., Markert, B., 2001. Biomonitoring d evaluation and assessment of heavy metal concentrations from two german moss monitoring surveys. In: Burga, C.A., Kratochwil, A. (Eds.), Biomonitoring: General and Applied Aspects on Regional and Global Scales, Tasks for Vegetation Science. Springer Netherlands, pp. 73e95.

pez, V., Bert, V., Dimitriou, I., Friesl-Hanl, W., Kidd, P., Mench, M., Alvarez-L o Herzig, R., Janssen, J.O., Kolbas, A., Müller, I., Neu, S., Renella, G., Ruttens, A., Vangronsveld, J., Puschenreiter, M., 2015. Agronomic practices for improving gentle remediation of trace element-contaminated soils. Int. J. Phytoremediation 17, 1005e1037. http://dx.doi.org/10.1080/15226514.2014.1003788. Leblond, S., Gombert, S., Colin, J.L., Losno, R., Traubenberg, C.R.D., 2004. Biological and temporal variations of trace element concentrations in the moss species Scleropodium purum (Hedw.) Limpr. J. Atmospheric Chem. 49, 95e110. http://

dx.doi.org/10.1007/s10874-004-1217-8. s, L., Lequy, E., Sauvage, S., Laffray, X., Gombert-Courvoisier, S., Pascaud, A., Galsomie Leblond, S., 2016. Assessment of the uncertainty of trace metal and nitrogen concentrations in mosses due to sampling, sample preparation and chemical analysis based on the French contribution to ICP-Vegetation. Ecol. Indic. 71, 20e31. http://dx.doi.org/10.1016/j.ecolind.2016.06.046. Markert, B.A., Breure, A.M., Zechmeister, H.G., 2003. Bioindicators and Biomonitors. Nagajyoti, P.C., Lee, K.D., Sreekanth, T.V.M., 2010. Heavy metals, occurrence and toxicity for plants: a review. Environ. Chem. Lett. 8, 199e216. http://dx.doi.org/ 10.1007/s10311-010-0297-8. €der, W., Steinnes, E., Uggerud, H.T., 2014. Nickel, S., Hertel, A., Pesch, R., Schro Modelling and mapping spatio-temporal trends of heavy metal accumulation in moss and natural surface soil monitored 1990e2010 throughout Norway by multivariate generalized linear models and geostatistics. Atmos. Environ. 99, 85e93. http://dx.doi.org/10.1016/j.atmosenv.2014.09.059. €rv, V., 1955. On the colloidal nature of peat-forming mosses. Arch. Soc. Zool. Puustja Bot. Fenn. Vanamo 9, 257e272. R Core Team, 2015. R:A language and environment for statistical computing. R Foundation for Statistical Computing. Vienna, Austria http://www.R-project. org.
Angel

ndez, J., Ramo
n Aboal, J., Carballeira, A., 2011. Study of the air Ares, A., Ferna quality in industrial areas of Santa Cruz de Tenerife (Spain) by active biomonitoring with Pseudoscleropodium purum. Ecotoxicol. Environ. Saf. 74, 533e541. http://dx.doi.org/10.1016/j.ecoenv.2010.08.019. Ruhling, A., Tyler, G., 1968. An ecological approach to lead problem. Bot. Not. 121, 21. 
, J., Suchara, I., Valiulis, D., 2009. Changes in total Sakalys, J., Kvietkus, K., Sucharova concentrations and assessed background concentrations of heavy metals in moss in Lithuania and the Czech Republic between 1995 and 2005. Chemosphere 76, 91e97. http://dx.doi.org/10.1016/j.chemosphere.2009.02.009. Sandilyan, S., Kathiresan, K., 2014. Decline of mangroves e a threat of heavy metal poisoning in Asia. Ocean. Coast. Manag. 102, 161e168. http://dx.doi.org/10.1016/ j.ocecoaman.2014.09.025. Part A. Sauvage, S., Plaisance, H., Locoge, N., Wroblewski, A., Coddeville, P., Galloo, J.C., 2009. Long term measurement and source apportionment of non-methane hydrocarbons in three French rural areas. Atmos. Environ. 43, 2430e2441. http://dx.doi.org/10.1016/j.atmosenv.2009.02.001. €der, W., Holy, M., Pesch, R., Harmens, H., Ilyin, I., Steinnes, E., Alber, R., Schro Aleksiayenak, Y., Blum, O., Cos¸kun, M., Dam, M., Temmerman, L.D., Frolova, M.,
ska, K., Jeran, Z., Korzekwa, S., Krmar, M., Frontasyeva, M., Miqueo, L.G., Grodzin  kovsk
Kubin, E., Kvietkus, K., Leblond, S., Liiv, S., Magnússon, S., Man a, B., €ni, L., Piispanen, J., Rühling, Å., Santamaria, J., Spiric, Z., Suchara, I., Tho Urumov, V., Yurukova, L., Zechmeister, H.G., 2010. Are cadmium, lead and mercury concentrations in mosses across Europe primarily determined by atmospheric deposition of these metals? J. Soils Sediments 10, 1572e1584. http:// dx.doi.org/10.1007/s11368-010-0254-y. €der, W., Pesch, R., Scho €nrock, S., Harmens, H., Mills, G., Fagerli, H., 2014. Schro Mapping correlations between nitrogen concentrations in atmospheric deposition and mosses for natural landscapes in Europe. Ecol. Indic. 36, 563e571. http://dx.doi.org/10.1016/j.ecolind.2013.09.013. €kinen, A., 1997. Reference materials for largeSteinnes, E., Rühling, Å., Lippo, H., Ma scale metal deposition surveys. Accreditation Qual. Assur 2, 243e249. http:// dx.doi.org/10.1007/s007690050141.
ndez, J.A., 2014. Use of a moss Varela, Z., Aboal, J.R., Carballeira, A., Real, C., Ferna biomonitoring method to compile emission inventories for small-scale industries. J. Hazard. Mater 275, 72e78. http://dx.doi.org/10.1016/ j.jhazmat.2014.04.061. WHO, 2014. Air Quality Deteriorating in Many of the World's Cities [WWW Document]. WHO. URL. http://www.who.int/mediacentre/news/releases/2014/ air-quality/en/. accessed 9.22.15. Wood, S., 2006. In: Chapman, R. (Ed.), Generalized Additive Models: an Introduction. Hall/CRC. Zuur, A.F., Ieno, E.N., Walker, N., Saveliev, A.A., Smith, G.M., 2009. Mixed Effects Models and Extensions in Ecology with R, Statistics for Biology and Health. Springer New York, New York, NY.