# Are Grimmia Mosses Good Biomonitors for Urban Atmospheric Metallic Pollution? Preliminary Evidence from a French Case Study on Cadmium

Bartholomé Vieille, Isabelle Albert, Sébastien Leblond, Florian Couvidat, Éric Parent, Caroline Meyer

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Bartholomé Vieille, Isabelle Albert, Sébastien Leblond, Florian Couvidat, Éric Parent, et al. Are Grimmia Mosses Good Biomonitors for Urban Atmospheric Metallic Pollution? Preliminary Evidence from a French Case Study on Cadmium. Atmosphere, 2021, 12, pp.491. 10.3390/atmos12040491. hal-03202170

HAL Id: hal-03202170

https://hal.science/hal-03202170

Submitted on 20 Sep 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Are *Grimmia* Mosses Good Biomonitors for Urban Atmospheric Metallic Pollution? Preliminary Evidence from a French Case Study on Cadmium

Bartholomé Vieille <fnref n="1" />,†, Isabelle Albert <fnref n="1" />,*,†, Sébastien Leblond <fnref n="2" />,†, Florian Couvidat <fnref n="3" /> and Éric Parent <fnref n="1" />
and Caroline Meyer <fnref n="2" />,*

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Vieille, B.; Albert, I.; Leblond, S.; Couvidat, F.; Parent, É.; Meyer, C. Are *Grimmia* Mosses Good Biomonitors for Urban Atmospheric Metallic Pollution? Preliminary Evidence from a French Case Study on Cadmium. *Atmosphere* 2021, 12, 491. https://doi.org/10.3390/atmos12040491

Academic Editors: Olivier Delhomme
and Lena Novack

Received: 5 March 2021

Accepted: 12 April 2021

Published: 14 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** Assessment of human exposure to atmospheric metals is a challenge, and mosses seem to be good biomonitors to help this purpose. Lacking roots, they are easy to collect and analyze. However, to our knowledge, no formal comparison was made between cadmium (Cd) measurements in *Grimmia* mosses and alternative forecasts of atmospheric Cd pollution as those produced by the CHIMERE chemistry transport model. This work aims at studying this link to improve further biomonitoring. We compare 128 Cd measurements in the cemetery mosses of Paris and Lyon metropolitan areas (France) to CHIMERE Cd atmospheric forecasts. The area to consider around the cemetery for the CHIMERE forecasts has been defined by Kendall rank correlations between both information sources—Cd in mosses and CHIMERE Cd forecasts—from different area sizes. Then, we fit linear models to those two data sets including step-by-step different sources of uncertainty. Finally, we calculate moss predictions to compare predictions and measurements in the two cities. The results show an apparent link between the Cd concentrations in mosses and CHIMERE Cd forecasts including in addition the same unique covariate, the moss support (grave or wall), in the two cities. However, this model cannot be directly transposed from region to region because the strength of the link appears to be regional.

**Keywords:** atmospheric metal in France; biomonitor; pollution; cadmium; moss; CHIMERE

## 1. Introduction

Since the 1960s, terrestrial mosses have been used for monitoring studies of atmospheric metallic deposition at different scales [1,2]. Their morphological and physiological properties make them excellent sensors of atmospheric contaminants in particulate or soluble form. The lack of a root and vascular system makes them dependent on the atmospheric deposition for their nutrients. Mosses efficiently uptake and retain their essential elements, but simultaneously capture many other elements such as metals [3,4]. Moss biomonitoring is a comparative method which makes it possible to rank the sites relative to one another according to the concentration value measured in the samples [5]. It is traditionally accepted that these values are quantitative estimates of atmospheric deposition. At the European scale, some studies, using Spearman rank correlations and multivariate decision trees, showed that atmospheric deposition modeled by the European Monitoring and Evaluation Programme (EMEP) is the main factor determining the accumulation of cadmium (Cd) and lead (Pb) in pleurocarpous mosses (*Pleurozium schreberi*, *Hylocomium splendens*, *Hypnum cupressiforme*, and *Pseudoscleropodium purum*) in rural areas and far from emission sources [6–8]. Similarly, in rural areas at the scale of certain countries (Germany,

<footnote n="1">¹ MIA-Paris-Mathématiques et Informatique Appliquées, INRAE-AgroParisTech-Université Paris-Saclay, CEDEX 05, 75231 Paris, France; bartholv@gmail.com (B.V.); eric.parent@agroparistech.fr (E.P.)</footnote>

<footnote n="2">² UMS 2006 Patrimoine Naturel, OFB/MNHN/CNRS, 12 Rue Buffon, 75005 Paris, France; sebastien.leblond@mnhn.fr</footnote>

<footnote n="3">³ INERIS, 60550 Verneuil-en-Halatte, France; florian.couvidat@ineris.fr</footnote>

<footnote>* Correspondence: isabelle.albert@inrae.fr (I.A.); caroline.meyer@mnhn.fr (C.M.)</footnote>

<footnote>† These authors contributed equally to this work.</footnote>
Switzerland, Bulgaria, and Spain, for example), regression or correlation approaches, linking the concentrations measured in pleurocarpous mosses to the deposition data modeled by EMEP for Cd and Pb, also demonstrated a moderate or high significant correlation between these matrices and that the value of this link depends on the country [9–12]. However, it is necessary to take these results with caution because the relationship between metal concentrations accumulated by mosses and atmospheric deposition is complex and dependent on several factors such as age of mosses or part of analyzed mosses [5,13].

Nevertheless, these previous studies on the main factors determining the accumulation of Cd in pleurocarpous mosses have been carried out in rural and forest areas and far from emissions sources. Moss data can be used in human health studies because they make it possible to overcome the limitations of data sets concerning atmospheric metals. According to the European Directive 2004/107/EC and 2008/50/EC, only five metals are statutorily followed by physicochemical methods, and the modeled deposit maps over a large territory therefore only concern a few elements [14,15]. Human cohort studies are mainly carried out in urban and peri-urban areas [16], and it is therefore necessary to adapt the protocols used in rural areas to urban areas. Few studies using mosses as a passive bioaccumulator (mosses growing on a given site are directly collected and analyzed [17]) have been carried out in urban areas, but most of them used moss species of the *Grimmia* Hedw. genus [18–20]. These species are abundant in urban and peri-urban environments on stone or cement-based supports and can live in contaminated environments.

In this study, we have more particularly studied the accumulation of Cd by species of the *Grimmia* genus in the urban area associated with two French metropolitan areas: Paris and Lyon. Indeed, Cd is an element of essentially anthropogenic origin emitted mainly by road traffic, industry, and the residential sector [21] and recognized as carcinogenic element [22]. Moreover, in rural areas, a recent study demonstrated a positive association between Cd and mortality [14]. Mosses are collected in cemeteries which represent relatively homogeneous areas and in which, depending on the management mode, mosses can remain for long periods [19]. Moreover, cemeteries share common traits in all cities although their characteristics can vary at the local level (proximity to a road, an industry, a crop, use of phytosanitary products, etc.). Natali et al. [19] demonstrated that in Île-de-France (Paris metropolitan area), *Grimmia pulvinata* (Hedw.) Sm. collected from cemeteries was an effective bioaccumulator for studying the spatial variability of metal deposits.

For the first time, we will study more precisely the relationship between the Cd accumulated by mosses of the genus *Grimmia* in urban area and atmospheric Cd concentrations modeled by CHIMERE [23]. CHIMERE is a chemistry transport model designed for regional atmospheric composition [24]. Pollutant emissions, advection and mixing, chemistry, and deposition are the key processes affecting the chemical concentrations outputs provided by the model. Many research applications (transport, mixing, turbulence, chemistry, etc.) as well as operational forecasts (air quality networks, etc.) are based on CHIMERE, working either online or offline. In this application, CHIMERE temporal and spatial grids were used to evaluate air concentration of Cd in the study areas of Paris and Lyon. Then, we will design and test statistical models explaining the relationship between atmospheric modeled Cd and accumulated Cd in the two biggest metropolitan areas of France: Paris and Lyon. These models will consider stationary variables such as, for example, the altitude of the site, the type of support on which the mosses were taken, the distance to the center of the zone, as well as the measurement of Cd uncertainties. Finally, we will study the best model defined for Paris metropolitan area on the metropolitan area of Lyon (and reversing the cities) to verify the potential regionalization of the responses of mosses to Cd metallic contaminant.
## 2. Materials and Methods

### 2.1. Cd Measurements in Mosses in Cemeteries of Paris and Lyon Metropolitan Areas

#### 2.1.1. Studied Sites

The studied areas (Figure 1) are defined according to the population density based on INSEE criteria. Briefly, 4 categories of spatial units are defined: residential zones (between 2000 and 5000 inhabitants), industrial or commercial activity zones (sparsely populated or with 1000 employees), sparsely populated recreation zones, and municipalities with have not been divided according to the previous criteria (https://www.insee.fr/fr/information/2383389 (accessed on 13 April 2021)). For these municipalities, their population density made it possible to divide them into 4 classes of high population density with very low population density. Each of the studied areas correspond to a circular surface whose center coincides with the barycentre of the continuous area of high and moderate population density (urban and peri-urban area). For Paris metropolitan area, the studied area is delimited by a circle with a radius of 49 km, the origin of which is located in Paris (Montparnasse Cemetery). This area covers all the departments of the so-called “Île de France” area. For Lyon metropolitan area, the studied area is delimited by a circle with a radius of 30 km, centered in Lyon city.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** The studied metropolitan areas: cemeteries sampled in the Paris and Lyon metropolitan areas (above, the municipality of Paris; below, the municipality of Lyon). Red point: mosses collected from graves; black point: mosses collected from walls.</figcaption></figure>

#### 2.1.2. Moss Sampling and Analysis

All moss samples in Paris and Lyon metropolitan area were, respectively, collected during 3 weeks in May 2018 and 15 consecutive days in June 2018. The moss species collected is *Grimmia pulvinata*. All collected samples are exclusively bryophytic colonies developed on substrates whose composition is cement, mortar, or concrete. Two types of support were sampled: tombs and the flat top of the perimeter walls. For a given site, all the moss colonies are taken from the same support (graves or walls). The walls were chosen as the preferential support for harvesting mosses. To avoid risks of elemental transfers by water runoff, attention was paid that no overhanging tree canopy or surface was present above sampled moss cushions. Mosses directly under the influence of metallic objects (flowerpots, religious objects, etc.) were rejected. In the Paris metropolitan area, mosses
were collected from graves in 58 cemeteries and from walls in 19 cemeteries. In the Lyon metropolitan area, mosses were collected from graves in 9 cemeteries and from walls in 42 cemeteries.

Each moss sample is made up of a minimum of 50 distinct colonies. The colonies are stored in a double-sealed zip bag and are stored in a cooler while awaiting return to the laboratory and placing them in a cold room at 4 °C. In the laboratory, they are dried at room temperature. Then, non-moss organic matter debris and gravels are removed. Whole shoots of moss samples are then ground with an automatic non-polluting titanium grinder (Pulverisette 14, Fritsch, Germany). To avoid any accidental contamination of our samples, they are taken and then handled with nitrile gloves and ceramic knives and a scalpel. A total of 250 mg of homogeneous powder is mineralized with a mixture of HF/HNO₃/H₂O₂. Cadmium (Cd) in mosses has been analyzed by inductively coupled plasma mass spectrometry (ICP-MS) by the USRAVE (INRAE- Centre de Bordeaux) (Limit of quantification is 0.02 µg/g). The device used by the USRAVE is an Agilent 7700x spectrometer. The ICP-MS is calibrated by external calibration, using multi-element solutions produced every 6 months. The internal calibration procedure involves the addition of the standard to obtain the same concentration in each sample, thus allowing the monitoring of the good performance of the device and the correction of drift over time. In this study, certified reference materials for trace elements (CRMs), namely, NIST 1547 (peach leaves), NIST 1573a (Tomato leaves), as well as European moss survey program reference materials—Moss M2 and Moss M3—were processed with each batch of samples. The Moss M2 and Moss M3 standards are powders of *Pleurozium schreberi* from two Finnish sites [25].

## 2.2. Atmospheric Cd Forecasts from Chimere Chemistry Transport Model

CHIMERE is an open source multi-scale mesoscale physico-chemistry model, jointly developed as a National Tool of the French Institut des Sciences de l'Univers, by IPSL/LMD (Palaiseau), INERIS (Verneuil en Halatte) and IPSL/LISA (Creteil) in France. The source code and the associated documentation are available on the website http://www.lmd polytechnique.fr/chimere (accessed on 13 April 2021). It simulates the troposphere (the first atmospheric layer of 10 km of altitude) with various horizontal resolutions over study areas ranging from city to continent. The CHIMERE model [24,26] solves integro-differential partial equations, coupled with meteorological forcing covariates and various sources of anthropogenic or natural emissions (Figure 2). The numerical code is parallelized and can deal with long-term simulations of the future behavior of many chemical elements, gases, and aerosols [27], as well as short-term operational forecasting of atmospheric pollution [28].

As main outputs, CHIMERE can deliver hourly pollutant concentrations (gases and aerosols) in 3D on any simulation domain. The various pollutants commonly monitored for atmospheric pollution can be simulated: O₃, NO, NO₂, SO₂, PM₂.₅, PM₁₀, and much more in dedicated CHIMERE versions, like some heavy metals (Pb, Cd, As, Ni, Cu, Zn, Cr, and Se) [26]. The retrospective work hereafter focuses on Cd: CHIMERE was therefore launched to provide monthly series of Cd concentrations in the air during year 2015 under the form of twelve 240 × 324 pixel grids covering the mainland of France in Europe.

In France, there are 4 main sources of Cd emissions: the manufacturing industry, energy transformation (incineration plants), road transport, and residential sector [21]. The CHIMERE model is based on the real emissions of Cd by each sector and industry, so each Cd data synthesizes the influence of each Cd source at a given pixel.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** General principle of a chemistry transport model such as CHIMERE. In the box "Meteorology", u* stands for the friction velocity, Q₀ is the surface sensible heat flux, and BLH is the boundary layer height. [c]ₘ₉₉ and [c]ₐ₉₉ are the modeled and the observed chemical concentrations fields, respectively. This figure is Figure 1.1 from https://www.lmd.polytechnique.fr/chimere/docs/CHIMEREdoc_v2020r1.pdf (accessed on 13 April 2021).</figcaption></figure>

### 2.3. Defining the Area of Influence of the Atmospheric Cd on Mosses

In a first exploratory part, we wanted to explore the distance from the cemetery to consider for establishing the link between the CHIMERE Cd atmospheric forecasts and the Cd measurements in mosses. We had one CHIMERE forecast for a squared zone of 3 km × 3 km, not a forecast at each cemetery site. Moreover, we did not know a priori the influence of reasonably far atmospheric metals on a cemetery site. Therefore, from each cemetery, we established different circular zones with a d radius, with the d distance varying from d = 2.5 to 50 km, and in each zone we calculated the mean of all the CHIMERE forecasts (one unique forecast in the area with d = 2.5 km). Then, for each d area, considering all cemeteries, we calculated the Kendall rank correlation between the measurements in mosses and the CHIMERE means attached to the d area around each cemetery. Finally, we had a Kendall rank correlation by d distance, and we could determine the distance that maximized the Kendall rank correlations between the CHIMERE forecasts and the measurements in mosses. The R (https://www.r-project.org/ (accessed on 13 April 2021)) package Kendall allows to calculate these Kendall rank correlations.

### 2.4. Linear Models with Fixed Effects

When observing a linear link between two variables of interest (here Cd measurement in mosses and CHIMERE Cd forecast), it is natural to produce linear models with other known explanatory covariates to explore in depth the relation between both variables. Here, we chose to explain Cd in mosses by the CHIMERE forecast with other covariates collected at the time of the moss sampling, such as the altitude of the sampling site, A; the type of support on which the moss was taken, S; the presence or not of traces of a phytosanitary treatment, T; or the transect of sampling to which the moss belongs, r.

Starting model: in this model, i indexes the mosses (i ∈ {1, ..., 77}, for example, in Paris region). The i's are assumed independent and identically distributed. The starting model is the following:

$$Cd_i = \alpha_{r(i)} + \beta_{r(i)}d_i + \gamma A_i + \zeta S_i + \chi C_i + \psi T_i + \epsilon_i,$$

where Cd is the Cd in the moss, d is the distance from the data sampling center (as explained previously the cemeteries were chosen from a central site and the observed Cd concentrations seem decreasing with the distance from this center), and Cᵢ is the CHIMERE Cd forecast mean from the 4 meshes around the cemetery i. The αᵣ(i) and βᵣ(i) are here
fixed effects depending on the sampling transects (8 sampling transects, see Figure 1). In addition, $\forall i, \text{Var}(\epsilon_i) = \sigma^2$, where $\sigma^2$ denotes the measure of uncertainty (variance; the same for each i moss).

Model selection: starting from this fairly rich model, one then proceeds to a selection of sub-models using the AIC criterion and a forward/backward method testing the sub-models [29]. The R function *lm* allows one to estimate the models, calculate the AIC criterion, and perform the forward/backward method.

## 2.5. Linear Models with Random Effects

As an extension of the linear model with fixed effects, the random effect linear model introduces additional random sources of variability within the data [30,31]. In addition to moss measurement errors, there may occur some randomness stemming from the sampling of transects. These reported transects are only a sample of transects chosen among a population of the many other possible ones: variability between transects has to be introduced into the statistical analysis. Because of this between-transect diversity, correlation (due to similarities between mosses records at sites belonging to the same transect) is expected when grouping measurements by transect. Then, we can no longer suppose that moss data are independent as thought previously (in the linear model with fixed effects), but they are now independent only conditionally to the knowledge of the effect due to their transect belonging.

Starting model: in this model, *i* indexes the mosses (*i* ∈ {1, ..., 77} for Paris region and *i* ∈ {1, ..., 51} for Lyon region). The *i*'s are independent and are assumed to be identically distributed conditionally to their transect. The starting model is the following:

$$Cd_{r(i)} = \alpha_{r(i)} + \beta_{r(i)}d_i + \gamma A_i + \zeta S_i + \chi C_i + \psi T_i + \epsilon_i,$$

where Cd is the Cd in the moss and the covariates are the same than in the fixed effects linear model but here the $\alpha_{r(i)}$ and $\beta_{r(i)}$ are random effects depending on the sampling transects (8 sampling transects as seen previously). Then,

$$\forall i, \text{Var}(\epsilon_i) = \sigma^2 + \sigma_r^2,$$

where $\sigma^2$ denotes the measure of uncertainty of the mosses and $\sigma_r^2$ denotes the measure of uncertainty of the transects.

Model selection: starting from this fairly rich model, one then proceeds to a selection of sub-models using the AIC criterion and a forward/backward method testing the sub-models [29]. The R function *lmer* from package *lme4* allows us to estimate the models and calculate the AIC criteria for the forward/backward method.

## 2.6. Linear Models with Experimental Errors

We also had an additional information in the mosses data set: the experimental measurement error was provided by the laboratory that made the quantification of Cd in mosses. Two standard deviations of this error amount to 20% of the Cd value measured in each moss. Thus, the greater the measured value, the greater the uncertainty. We will therefore take into account this measurement error which can have a strong influence on the statistical fit.

Starting model: in this model, the *i*'s are independent but are no longer assumed to be identically distributed: they have a variance proportional to the measurement error provided. The starting model is the following:

$$Cd_i = \alpha_{r(i)} + \beta_{r(i)}d_i + \gamma A_i + \zeta S_i + \chi C_i + \psi T_i + \epsilon_i,$$
for $i \in \{1, \dots, 77\}$ (Paris) and $i \in \{1, \dots, 51\}$ (Lyon), where the same notations as previously are adopted. The $\alpha_{r(i)}$ and $\beta_{r(i)}$ are here fixed effects because previous results did not justify the choice of a random model for them. In addition,

$$ \forall i, \mathbf{Var}(\epsilon_i) = \sigma^2 \frac{\xi_i^2}{4}, $$

where $\xi_i$ denotes the measurement error provided by the laboratory for the $i$-th moss and $\sigma^2$ denotes the residual measure of uncertainty of the mosses.

Model selection: starting from this fairly rich model, one then proceeds to a selection of sub-models using the AIC criterion and a forward/backward method testing the sub-models [29]. The R package gls allows to estimate the models and perform the forward/backward method testing the sub-models.

## 2.7. Linear Model with Chimere Errors

So far, we had considered the variable CHIMERE $C_i$ as a deterministic explanatory quantity. However, remember that these data come from an interpolation of the meshes whose centers are close to the moss site. Indeed, we did not have CHIMERE data located at exactly the same points as the mosses but only a network of predictions that we had assumed to be uniform on each cell. In the following, we propose to interpolate the CHIMERE data at mosses points in another way, using two mathematical tools which are the variogram and the kriging [32]. We therefore need to introduce a new uncertainty linked to the kriging estimation error.

As the previously introduced model selects the same variables whether for Lyon or Paris, we did not make a selection of sub-models but we relied on this same choice of covariates. Additionally, we will suppose that the Cd concentration CHIMERE $C_i$ is itself a random process following at each point of Île-de-France and Lyon a normal law. We get the expectation (noted $C_i^*$) and the variance (noted $V_i^2$) of this process in each of the moss points ($i \in \{1, \dots, 77\}$ for Paris) thanks to the kriging carried out at the moss points using the CHIMERE data. The regression problem is now written mathematically as follows:

$$ \left\{ \begin{array}{l} \mathrm{Cd}_i = \alpha + \chi C_i + \zeta S_i + \epsilon_i, \\ \epsilon_i | C_i \sim \mathcal{N}(0, \sigma^2 \xi_i^2), \\ C_i \sim \mathcal{N}(C_i^*, V_i^2), \end{array} \right. $$

where $\sigma^2$ is a value to estimate which allows a maximum of flexibility to the error uncertainty of each moss. $\xi_i$ is still the measurement error provided by the laboratory for the $i$-th moss.

**Kriging at moss points:** We started by constructing a variogram on the CHIMERE data. The empirical variogram was then fitted using a Matérn variogram model [32]. We then used kriging at the blunt points (77 for Paris and 51 for Lyon). This kriging gave us at each of these points an average and a variance which therefore correspond to $C_i^*$ and $V_i$, respectively.

**Parameter estimation:** As there was no predefined function in R to deal with this type of problem, we chose to solve the previous optimization problem “manually” by maximizing the log-likelihood. According to the model stated just above, $Cd_i$ marginally follows a normal distribution, $N(\alpha + \chi C_i^* + \zeta S_i, \chi^2 V_i^2 + \sigma^2 \xi_i^2)$. Moreover, the $Cd_i$ will be assumed to be independent, even if there is a small covariance between the $C_i$. Therefore, the log-likelihood of the n-sample is given by

$$ \ln(L(\mathrm{Cd}_1, \dots, \mathrm{Cd}_n)) = \sum_{i=1}^{n} \ln\left(\frac{1}{\sqrt{(2\pi)(\chi^2 V_i^2 + \sigma^2 \xi_i^2))}}\right) - \frac{1}{2} \frac{(\mathrm{Cd}_i - (\alpha + \chi C_i^* + \zeta S_i))^2}{\chi^2 V_i^2 + \sigma^2 \xi_i^2}. $$

Dropping the constant $2\pi$ term, this leads to
$$ \ln(L(\mathrm{Cd}_1, \dots, \mathrm{Cd}_n)) = -\frac{1}{2} \left(\sum_{i=1}^{n} \ln(\chi^2 V_i^2 + \sigma^2 \xi_i^2) + \frac{\mathrm{Cd}_i - (\alpha + \chi C_i^* + \zeta S_i)^2}{\chi^2 V_i^2 + \sigma^2 \xi_i^2}\right). $$

We then maximized this log-likelihood using the BFGS algorithm with the R package *lbfgs*. The value returned by this function corresponds to $-2 \times \ln(L)$, where L denotes here the likelihood. The parameters were obtained at the end of this optimization. Their confidence intervals were obtained by bootstrap by fitting the model on 6000 samples obtained by drawing with replacement on data from mosses from the Paris and Lyon region.

## 2.8. Predictions at One City from a Model Fitted to the Other One

Finally, we wanted to explore the predictive value of the best founded model regarding the AIC criterion. For this, we calculated the moss measurements predictions at Lyon cemeteries from the best model fitted at Paris (we use the values of the Lyon covariates, the estimated values of the parameters from the model fitted at Paris and the model equation). Then, we compare those predictions with the observed measurements at Lyon to explore the predictive value of the best model. Furthermore, we did the same thing reversing the cities.

# 3. Results

## 3.1. Characterization of Distribution of Cd Accumulated by Mosses and Cd Modeled by Chimere

The mean value of the accumulated Cd in *Grimmia* is higher for the Paris metropolitan area (0.99 ± 0.56 µg/g for a range between 0.25 µg/g and 2.79 µg/g) than for Lyon metropolitan area (0.41 ± 0.26 µg/g for a range between 0.15 µg/g and 1.27 µg/g) (Table 1). However, the Lyon metropolitan area has a greater variability in Cd concentration than the Paris metropolitan area (64% and 57%, respectively) as well as a greater number of atypical values. The high atypical threshold calculated in the Paris metropolitan area is 2.48 µg/g, and two sites have a concentration value greater than this threshold (2.6% of sites). High concentrations of Cd are measured in cemeteries located in the inner suburbs of Paris as well as along the northern axis (Figure 3). The high atypical threshold calculated for Lyon metropolitan area is 0.80 µg/g and five sites have a concentration value greater than this threshold (9.8% of sites). The strong Cd values are located in the city of Lyon as well as in the east of the study area (Figure 3).

<table><caption>**Table 1.** Main statistics of Cd in mosses (µg/g).</caption><thead><tr><td>Parameters</td><td>Value for Paris</td><td>Value for Lyon</td></tr></thead><tbody><tr><td>Mean</td><td>0.99</td><td>0.41</td></tr><tr><td>Standard deviation</td><td>0.56</td><td>0.26</td></tr><tr><td>Min</td><td>0.25</td><td>0.15</td></tr><tr><td>Max</td><td>2.79</td><td>1.27</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Maps showing the value of Cd (µg/g) in mosses. On the left, Paris region; on the right, region of Lyon.</figcaption></figure>
Figure 4 shows the distribution of Cd modeled by CHIMERE in the meshes of the Paris and Lyon regions. Table 2 gives the main descriptive statistics of these values in the two regions.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Maps showing the value of the logarithm of the concentration in air of CHIMERE Cd (log(µg/m³)). On the **left**, Paris region; on the **right**, region of Lyon.</figcaption></figure>

<table><caption>**Table 2.** Main statistics of Cd concentrations (µg/m³) modeled by CHIMERE.</caption><thead><tr><th>Parameters</th><th>Value for Paris</th><th>Value for Lyon</th></tr></thead><tbody><tr><td>Mean</td><td>3.55 × 10<sup>-5</sup></td><td>3.47 × 10<sup>-5</sup></td></tr><tr><td>Standard deviation</td><td>1.64 × 10<sup>-5</sup></td><td>1.41 × 10<sup>-5</sup></td></tr><tr><td>Min</td><td>2.24 × 10<sup>-5</sup></td><td>1.79 × 10<sup>-5</sup></td></tr><tr><td>Max</td><td>2.61 × 10<sup>-4</sup></td><td>1.24 × 10<sup>-4</sup></td></tr></tbody></table>

From the maps of this section, one observes a link between Cd in mosses and CHIMERE Cd. In fact, whether for the region of Lyon or the region of Paris, the Cd in mosses is highest where the concentration in the air of CHIMERE Cd is the highest: whether for mosses or CHIMERE meshes, Cd is much more present in the heart of the two regions.

## 3.2. The Area of Influence of the Atmospheric Cd on Mosses

We calculated the Kendall rank correlation between the Cd in mosses and the Cd concentration in the atmosphere averaged in a circle of a *d* radius around cemeteries: as a function of the distance *d* for Paris and Lyon metropolitan areas, the values of the Kendall rank correlations showed differences between areas (Figure 5). For Paris metropolitan area, the Kendall rank correlation is maximal (0.55) when *d* is below 5 km and minimal (0.44) when *d* is equal 50 km (Paris metropolitan area was a circle with a radius of 49 km when defining moss sampling). However, the Kendall rank correlation level is very stable before 40 km (Figure 5). For the Lyon metropolitan area, the Kendall rank correlation is maximal (0.49) when *d* is equal to 25 km and minimal (0.18) when *d* is equal to 50 km (Figure 5). The decrease in the value of the correlation from 30 km can be explained by the fact that the size of the study area for moss sampling is 30 km for Lyon metropolitan area. In addition, for Lyon metropolitan area, Kendall rank correlation level is very stable below 30 km. There are some oscillations, but we must not forget that the level of correlation is already high and these oscillations only reveal small variations in correlations.

As the moss can be located at the very edge of the mesh, we have chosen in the following to consider the average of the concentration in the air on the four meshes whose centers are closest to the moss. The weighting of each of the considered meshes is inversely proportional to the distance which separates the moss from the core of the CHIMERE mesh.
<figure><img src="image_9.png" /><figcaption>Figure 5. Evolution of the Kendall rank correlation as a function of the distance from which the concentration in air of cadmium is averaged around the moss. On the left, the evolution in the Paris region, and on the right, the evolution in the Lyon region.</figcaption></figure>

Figure 6 represents the Cd of mosses as a function of the CHIMERE Cd averaged as proposed. We clearly identify a linear link (clearer in Paris) which suggests the use of the linear model in the following.

<figure><img src="image_10.png" /><figcaption>Figure 6. Link between Cd in mosses and the average CHIMERE Cd of the closest neighbors, for Île-de-France and the Lyon region.</figcaption></figure>

### 3.3. Linear Models with Fixed Effects

For the Paris data fitting, the model obtained after the model selection is as follows: for $i \in \{1, \dots, 77\}$:

$$Cd_i = \alpha_{r(i)} + \zeta S_i + \chi C_i + \gamma A_i + \epsilon_i,$$

where $Cd_i$ is the Cd in the moss $i$, $S_i$ is the type of support on which the moss was taken (grave or wall), $C_i$ is the CHIMERE Cd forecast mean from the four meshes around the moss $i$, and $A_i$ is the altitude of the sampling site. The $\alpha_{r(i)}$ are fixed intercepts depending on the sampling transects. In addition, $\forall i, \text{Var}(\epsilon_i) = \sigma^2$, where $\sigma^2$ denotes the measure of uncertainty (the same for each $i$-th moss). The choice of this model decreases the AIC by 10 compared to the initial model, whose AIC was 81. When selecting the model, we used the backward-forward method [29]; this means that all the sub-models, including those having only one explanatory variable or not containing the support, are less relevant.

For the Lyon data fitting, the model obtained after the model selection is as follows: for $i \in \{1, \dots, 51\}$:

$$Cd_i = \alpha + \zeta S_i + \chi C_i + \epsilon_i,$$
where Cd$_i$ is the Cd in the moss i, S$_i$ is the type of support on which the moss was taken (grave or wall), and C$_i$ is the CHIMERE Cd forecast mean around the moss i. The α is a fixed intercept not dependent on the sampling transects. In addition, ∀i, **Var**($\epsilon_i$) = σ², where σ² denotes the measure of uncertainty (the same for each i-th moss). The AIC decreases by 18 compared to the initial model, whose AIC was -6.

### 3.4. Linear Models with Random Effects

The model obtained after the model selection for Paris is as follows: for $i \in \{1, \dots, 77\}$:

$$Cd_i = \alpha_{r(i)} + \zeta S_i + \gamma A_i + \chi C_i + \epsilon_i,$$

where Cd$_i$ is the Cd in the moss i, S$_i$ is the type of support on which the moss was taken (grave or wall), C$_i$ is the CHIMERE Cd forecast mean around the moss i, and A$_i$ is the altitude of the sampling site. The $\alpha_{r(i)}$ are fixed intercepts depending on the sampling transects. In addition, ∀i, **Var**($\epsilon_i$) = σ², where σ² denotes the measure of uncertainty (the same for each i-th moss). The choice of this model decreases the AIC by 5.5 compared to the initial model, whose AIC was 76.5. The backward-forward method has also been used [29].

For Lyon, $i \in \{1, \dots, 51\}$:

$$Cd_i = \alpha + \zeta S_i + \chi C_i + \epsilon_i,$$

where Cd$_i$ is the Cd in the moss i, S$_i$ is the type of support on which the moss was taken (grave or wall), and C$_i$ is the CHIMERE Cd forecast mean around the moss i. The α is a fixed intercept not dependent on the sampling transects. In addition, ∀i, **Var**($\epsilon_i$) = σ², where σ² denotes the measure of uncertainty (the same for each i-th moss). The AIC criterion decreases by 5 from -19 to -24.

We observe that the selected model is the same for Lyon data and Paris data than these in the previous section: no longer containing random effects.

### 3.5. Linear Models with Experimental Errors

The model obtained after the model selection is as follows: for $i \in \{1, \dots, 77\}$ (for Paris) and $i \in \{1, \dots, 51\}$ (for Lyon):

$$Cd_i = \alpha + \zeta S_i + \chi C_i + \epsilon_i,$$

where Cd$_i$ is the Cd in the moss i, S$_i$ is the type of support on which the moss was taken (grave or wall), and C$_i$ is the CHIMERE Cd forecast mean around the moss i. The α is a fixed intercept not dependent on the sampling transects. In addition, ∀i, **Var**($\epsilon_i$) = σ²ζ²/4, where ζ$_i$ denotes the measurement error provided by the laboratory for the i-th moss and σ² denotes the residual measure of uncertainty. For the Paris data fitting, the choice of this model decreases the AIC by 231.8 compared to the initial model, whose AIC was 236.7. The backward-forward method has also been used [29]. For the Lyon data fitting, the same model is chosen and the AIC criterion decreases by 249, from 180 to -69. In both cases, the difference in AIC is very important.

### 3.6. Linear Models with Experimental Errors and Chimere Errors

Given the previous results, we fitted the last model selected in both cities adding only uncertainty on the CHIMERE forecasts (see Section 2.7 for the complete model formulation). For the Paris metropolitan area, the value of AIC criterion is -132.79 with 5 degrees of freedom. For the Lyon metropolitan area, the value of AIC is -160.42 with 5 degrees of freedom. The backward-forward method has also been used [29]. This model is the one with the smallest AIC among the models seen previously. It therefore appears to be our best model, from which we should make our predictions.
Table 3 presents the value of the parameters of the model, adjusted for Lyon and Paris. Note that the value of the coefficient “support” ($\zeta$) is positive indicating that the Cd in the moss is higher on the graves ($S_i = 1$).

### 3.7. Predictions

In this part, we present the predictions in the Lyon (resp. Paris) metropolitan area for the preferred model, namely, the model with experimental errors and CHIMERE errors, fitted on the Paris (resp. Lyon) metropolitan area (Figure 7). The two models have difficulties fitting the high values of Cd even if their large confidence intervals for higher Cd values almost encompass the observed values. Figure 8 shows the predictions in Lyon (on the left) and in Paris (on the right) if the model is fitted in the other city (Paris on the left and Lyon on the right) for evaluating the universality of the two models. Here, the predictions may appear poor for the small observed Cd values in Lyon because the model overfits those values, and for the observed Cd values in Paris, they are all in the intervals of confidence given by the Lyon's model but these intervals are very large. These results are not surprising looking at the value of the $\chi$ parameter in the model fitted to Paris compared to the value of the $\chi$ parameter in the model fitted to Lyon (Table 3).

<figure><img src="image_11.png" /></figure>

**Figure 7.** Predictions of Cd in mosses by the model taking into account experimental errors and CHIMERE errors, as well as their confidence intervals. On the left (resp. right): prediction of Cd in Paris (resp. Lyon) mosses by the model trained on Paris (resp. Lyon) mosses. The line corresponds to the line $y = x$. Red (resp. Blue) is used when the support of the corresponding moss is a wall (resp. a grave).

<figure><img src="image_12.png" /></figure>

**Figure 8.** Predictions of Cd in mosses by the model taking into account the measurement uncertainty provided as well as their confidence intervals. On the left (resp. right): prediction of cadmium in Lyon (resp. Paris) mosses by the model trained on Paris (resp. Lyon) mosses. The line corresponds to the line $y = x$. Red (resp. Blue) is used when the support of the corresponding moss is a wall (resp. a grave).
<table><caption>**Table 3.** Model parameter estimations: mean values and 95% confidence intervals obtained by bootstrap for Paris and Lyon metropolitan areas. The number of bootstrap samples is 6000.</caption><thead><tr><th>Parameters</th><th>α̂</th><th>λ̂</th><th>ζ̂</th><th>σ̂</th></tr></thead><tbody><tr><td>Mean value for Paris</td><td>0.02</td><td>9000</td><td>0.22</td><td>1.39</td></tr><tr><td>Mean value for Lyon</td><td>0.12</td><td>2778</td><td>0.11</td><td>1.65</td></tr><tr><td>Confidence interval (95%) for Paris</td><td>[-0.04, 0.11]</td><td>[8936, 9175]</td><td>[0.10, 0.31]</td><td>[1.10, 1.62]</td></tr><tr><td>Confidence interval (95%) for Lyon</td><td>[0.09, 0.15]</td><td>[2476, 3370]</td><td>[0.03, 0.33]</td><td>[1.35, 1.84]</td></tr></tbody></table>

## 4. Discussion

The objectives of this study were to better understand the relationship between the atmospheric Cd concentrations modeled by CHIMERE in two French urban areas and Cd concentrations accumulated by mosses of the genus *Grimmia* in the cemeteries of these two metropolitan areas. Although Paris is more densely populated than Lyon, the traffic is denser and the most emitting sites are located there, the two metropolises show similar average Cd concentrations (Table 2), and a higher dispersion coefficient is found in Paris urban area (46%) than in Lyon urban area (41%). This global result is not surprising because the emission sources between the two cities are the same (energy transformation, road transport, and residential sector [21]) but the size of studied area is different. Thus, for the Paris urban area, the study area includes more non-highly urbanized areas than the Lyon urban area. Nevertheless, the minimum and maximum values are higher in Paris than in Lyon. However, these values remain lower than those measured from PM in other European cities. For example, atmospheric Cd concentrations ranged from $5.3 \times 10^{-4}$ to $9.2 \times 10^{-4}$ µg/m³ in Warsaw in 2017 [33] and were significantly lower than those measured in southern China (0.019 µg/m³) [34]. The authors of two recent studies used *Grimmia* mosses as a passive bioaccumulator of metallic elements in the urban environment [19,20]. In our study, the concentrations of Cd accumulated by *Grimmia* mosses are equivalent to those found by Gallego-Cartagena et al. in north Spain [20] in sites under the influence of road traffic and industry, which could demonstrate that the origin of the contamination is atmospheric and not due to an effect of variables such as the substrate of mosses. The study carried out in Île-de-France did not analyze the Cd contents in mosses [19], but demonstrated that *Grimmia pulvinata* is a good proxy for atmospheric pollution by 20 elements (Al, As, Br, Ca, Ce, Cl, Cr, Cu, Fe, K, Mn, Ni, P, Pb, Rb, S, Sr, Ti, and Zn). This study highlights differences in the level of contamination in urban areas dependent on the distance from the urban environment [19]. Our mosses sampled in the cemeteries of the Paris and Lyon metropolitan areas have Cd contents which vary in the study area according to the atmospheric concentrations of Cd modeled by CHIMERE model. Our study makes it possible to add Cd to this list of elements studied in Île-de-France and makes it possible to specify that the variations of concentrations in the mosses are linked to the variations of concentrations in the atmosphere in two metropolitan areas of the territory.

Depending on the chosen spatial resolution, the Kendall correlations between the two series range between 0.44 and 0.56 for Paris metropolitan area and between 0.15 and 0.49 for Lyon metropolitan area. These rather weak correlations may stem from the fact that the available Cd atmospheric values simulated by the CHIMERE model date back to 2015, i.e., 3 years before the mosses were collected. However, the study of Cd emitted values since 1990 (from emission data compiled by European Monitoring and Evaluation Program (EMEP)) shows that they have remained stable since 2010. Moreover, to our knowledge, there is no other study on correlations between atmospheric Cd modeled on an urban area and the values of Cd concentrations in acrocarp mosses of the genus *Grimmia*. Cd concentrations in mosses have previously been correlated with values of atmospheric deposition using EMEP transport models [6,7]. In these studies, the correlation values are significantly higher. Furthermore, Berg et al. [35] and Berg and Steinnes [36] find correlations of $r^2 = 0.82 - 0.84$ between the Cd concentrations in *Hylocomium splendens* and *Pleurozium schreberi* and the wet deposits of Cd for mosses sampled as part of the European
monitoring program. In mosses sampled in Europe in 2005/2006, Schröder et al. [9] found a correlation of 0.65 between the Cd deposition of the sampling year and the concentrations in mosses and of 0.63 when considering the deposit over 3 years (Spearman correlations). This same study also showed that Cd concentrations in mosses were correlated with regional characteristics of the study area (urban land use at different radius).

Contrary to these studies, our results did not show a correlation between the Cd deposits modeled by CHIMERE near cemeteries and our Cd concentrations accumulated by *Grimmia* in cemeteries (results not shown). These results could be explained by the fact that the atmospheric deposition calculated by CHIMERE in a mesh 3 km × 3 km relies on the main land use in the mesh and this main land use, in an urban area, must be hugely different from the land use in the cemeteries. This difference of land uses should evidently induce a discrepancy between the modeled deposition in the meshes and what is measured in the cemeteries. In Germany, Nickel et al. [37] also demonstrated the importance of the deposition model taken into consideration in the correlations and showed a better correlation between modeled atmospheric deposition and Cd concentrations in mosses with deposit data modeled by LOTOS-EUROS which considers the land use, unlike deposits modeled by EMEP.

As for the studies carried out on mosses collected as part of the European monitor-
ing program [6,9,38], the values of Cd concentrations in mosses are correlated with the
concentrations calculated in the atmosphere near cemeteries but, as in these studies, the
atmospheric concentrations of Cd modeled cannot alone explain the Cd contents measured
in our samples of *Grimmia*. The values of Kendall correlations are lower in the metropolitan
area of Lyon than in that of Paris and suggest that regional factors are involved in the
relationship between atmospheric Cd and Cd accumulated by mosses [13].

The various regression models fitted in this study allowed us to explore this rela-
tionship in more details. First, before taking into account the uncertainty associated with
analytical measurement errors, the two types of adjusted models (fixed and random effects)
show that altitude is one of the factors influencing the results in the Paris metropolitan
area, as well as the sampling transects. However, once we consider the uncertainty related
to the analytical measurement (20% of the measurement), the model has recourse to the
same explanatory covariates for the two metropolitan areas, with AICs which suggest a
clear improvement in the fit of the models. In addition to Cd CHIMERE forecasts, only
the support on which moss was taken remains an explanatory variable in the model of the
two cities. Table 3 of the most complex model shows that the coefficient of this “support”
variable is not significantly different between metropolitan areas. In our study, the concen-
trations of Cd accumulated by the mosses collected from the graves are higher than for the
mosses collected from the walls of the cemetery. However, it is not possible to conclude
with certainty on the origin of this difference in concentrations between mosses taken from
graves and walls. Indeed, due to the religious character of the cemeteries and the number
of graves and walls sampled, it was impossible to take samples from the graves’ substrates
for analysis to crosscheck for diffusion. Nevertheless, a simple chemical analysis of the
cement would not have provided a definitive answer because the essential factor is the
bioavailability of the element in the substrate, and the metal content of the substrate is not
indicative of the values found in mosses [39]. This last model tested (linear model with
experimental errors and CHIMERE errors) based on a kriging, i.e., a linear combination
of all the surrounding CHIMERE cells whose coefficients are the correlations with the
cemetery, still allows a decrease in AIC and therefore even better overall adjustment in
Paris and Lyon metropolitan areas. This more complex consideration of the CHIMERE
data and their error makes it possible to take better account of the fact that these data
are imprecise at the mosses sampling coordinates because they are evaluated at another
geographical scale and the model is less accurate for high values of Cd (Figure 7). In the
Lyon metropolitan area, some average predictions are far apart from those measured, but
the intervals of uncertainties are important and the model is reliable. Yet, Table 3 shows
that the slope of the CHIMERE variable is significantly different in the two metropolitan
areas, thus reflecting a regional effect of the atmospheric Cd concentration. Therefore, Figure 8 shows that the two models are not reliable to be transposed to the other region. Further investigations should be performed to understand whether the introduction of regional effects could permit substantial improvements in fitting and in predictive use of such models.

**5. Conclusions**

Mosses of the genus *Grimmia* are efficient Cd bioaccumulators in urban and peri-urban environments and can be used as a proxy for atmospheric Cd concentrations. For the first time, our results show that the Cd concentrations in acrocarpous mosses are correlated with the modeled Cd concentrations by the CHIMERE chemistry transport model in two French metropolitan areas: Paris and Lyon. The better model is a linear model with experimental errors and CHIMERE errors based on a kriging and including the variable "substrate of moss sampling". The relationships between Cd concentrations in mosses and modeled Cd concentrations by the CHIMERE depend on some regional factors and the model developed for one metropolitan area is not fully reliable to predict the Cd observations of the other metropolitan area.

Despite the large difference in the coefficients of the model adjusted to Lyon or Paris, we obtain a consistent prediction of the level of cadmium present in the moss. However, in addition to the difference in coefficients between the two cities, it seems reasonable to think that this model can only be applied to large cities, where cadmium pollution is quite high, as is the case for Lyon and Paris. In rural areas, for example, initial investigations led us to believe that there was no such link on this scale. Studies in other urban and peri-urban areas would provide a better understanding of the regional variability of the models. Finally, it would be interesting to develop similar models for other elements modeled by CHIMERE and analyzed in our moss samples such as As, Ni, and Hg.

**Author Contributions:** Conceptualization, S.L. and C.M.; Data curation, F.C. and C.M.; Formal analysis, B.V., I.A., F.C. and É.P.; Funding acquisition, S.L. and C.M.; Investigation, B.V., É.P. and C.M.; Methodology, I.A., S.L., É.P. and C.M.; Project administration, S.L.; Supervision, I.A., F.C. and É.P.; Validation, I.A., S.L. and É.P.; Writing—original draft, B.V., I.A., É.P. and C.M.; Writing—review & editing, I.A., S.L., É.P. and C.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is a contribution of the ANR PoCoMo project (Air pollution effects on cognition in the CONSTANCES cohort). PoCoMo is supported by the French National Agency for Research under the "Vie, santé et bien-être" Program (ANR-17-CE36-0005). Further funding was provided to C.M. by the Fondation de France.

**Data Availability Statement:** The data of Cd accumulated by mosses presented in this study are available on request from the corresponding author (caroline.meyer@mnhn.fr).

**Acknowledgments:** The authors thank the municipalities, those in charge of the cemeteries, and the guards for their authorizations to take samples. The authors also thank Emeline Lequy-Flahault for its involvement and assistance in setting up the sampling strategy and Nora Rouillier for the preparation of moss samples.

**Conflicts of Interest:** The authors declare no conflict of interest.

**References**

1. Onianwa, P. Monitoring atmospheric metal pollution: A review of the use of mosses as indicators. *Environ. Monit. Assess.* **2001**, *71*, 13–50. [CrossRef] [PubMed]

2. Harmens, H.; Norris, D.; Steinnes, E.; Kubin, E.; Piispanen, J.; Alber, R.; Aleksiayenak, Y.; Blum, O.; Coşkun, M.; Dam, M.; et al. Mosses as biomonitors of atmospheric heavy metal deposition: Spatial patterns and temporal trends in Europe. *Environ. Pollut.* **2010**, *158*, 3144–3156. [CrossRef] [PubMed]

3. Bates, J. Mineral nutrient acquisition and retention by bryophytes. *J. Bryol.* **1992**, *17*, 223–240. [CrossRef]

4. Chakrabortty, S.; Paratkar, G.T. Biomonitoring of trace element air pollution using mosses. *Aerosol Air Qual. Res.* **2006**, *6*, 247–258. [CrossRef]
5. Aboal, J.; Fernández, J.; Boquete, T.; Carballeira, A. Is it possible to estimate atmospheric deposition of heavy metals by analysis of terrestrial mosses? *Sci. Total. Environ.* **2010**, *408*, 6291–6297. [CrossRef] [PubMed]

6. Holy, M.; Pesch, R.; Schröder, W.; Harmens, H.; Ilyin, I.; Alber, R.; Aleksiayenak, Y.; Blum, O.; Coşkun, M.; Dam, M.; et al. First thorough identification of factors associated with Cd, Hg and Pb concentrations in mosses sampled in the European Surveys 1990, 1995, 2000 and 2005. *J. Atmos. Chem.* **2009**, *63*, 109–124. [CrossRef]

7. Schröder, W.; Holy, M.; Pesch, R.; Harmens, H.; Ilyin, I.; Steinnes, E.; Alber, R.; Aleksiayenak, Y.; Blum, O.; Coşkun, M.; et al. Are cadmium, lead and mercury concentrations in mosses across Europe primarily determined by atmospheric deposition of these metals? *J. Soils Sedim.* **2010**, *10*, 1572–1584. [CrossRef]

8. Harmens, H.; Ilyin, I.; Mills, G.; Aboal, J.; Alber, R.; Blum, O.; Coşkun, M.; De Temmerman, L.; Fernández, J.; Figueira, R.; et al. Country-specific correlations across Europe between modelled atmospheric cadmium and lead deposition and concentrations in mosses. *Environ. Pollut.* **2012**, *166*, 1–9. [CrossRef] [PubMed]

9. Schröder, W.; Pesch, R. Long-term monitoring of the metal accumulation in forests measured by use of the moss technique. *Eur. J. For. Res.* **2010**, *129*, 475–488. [CrossRef]

10. Thöni, L.; Yurukova, L.; Bergamini, A.; Ilyin, I.; Matthaei, D. Temporal trends and spatial patterns of heavy metal concentrations in mosses in Bulgaria and Switzerland: 1990–2005. *Atmos. Environ.* **2011**, *45*, 1899–1912. [CrossRef]

11. Schröder, W.; Holy, M.; Pesch, R.; Zechmeister, H.; Harmens, H.; Ilyin, I. Mapping atmospheric depositions of cadmium and lead in Germany based on EMEP deposition data and the European Moss Survey 2005. *Environ. Sci. Eur.* **2011**, *23*, 19. [CrossRef]

12. Boquete, M.; Fernández, J.; Carballeira, A.; Aboal, J. Relationship between trace metal concentrations in the terrestrial moss Pseudoscleropodium purum and in bulk deposition. *Environ. Pollut.* **2015**, *201*, 1–9. [CrossRef] [PubMed]

13. Boquete, M.T.; Ares, A.; Fernández, J.; Aboal, J. Matching times: Trying to improve the correlation between heavy metal levels in mosses and bulk deposition. *Sci. Total. Environ.* **2020**, *715*, 136955. [CrossRef] [PubMed]

14. Lequy, E.; Siemiatycki, J.; Leblond, S.; Meyer, C.; Zhivin, S.; Vienneau, D.; de Hoogh, K.; Goldberg, M.; Zins, M.; Jacquemin, B. Long-term exposure to atmospheric metals assessed by mosses and mortality in France. *Environ. Int.* **2019**, *129*, 145–153. [CrossRef]

15. Meyer, C.; Leblond, S.; Jacquemin, B.; Lequy, E. Métaux, pollution de l'air et santé-Les mousses, des alliées originales en épidémiologie. *MÉDecine Sci.* **2020**, *36*, 376–381. [CrossRef]

16. Bentayeb, M.; Stempfelet, M.; Wagner, V.; Zins, M.; Bonenfant, S.; Songeur, C.; Sanchez, O.; Rosso, A.; Brulfert, G.; Rios, I.; et al. Retrospective modeling outdoor air pollution at a fine spatial scale in France, 1989–2008. *Atmos. Environ.* **2014**, *92*, 267–279. [CrossRef]

17. Markert, B. Definitions and principles for bioindication and biomonitoring of trace metals in the environment. *J. Trace Elem. Med. Biol.* **2007**, *21*, 77–82. [CrossRef] [PubMed]

18. Ugur, A.; Özden, B.; Saç, M.; Yener, G.; Altinbaş, Ü.; Kurucu, Y.; Bolca, M. Lichens and mosses for correlation between trace elements and 210Po in the areas near coal-fired power plant at Yatağan, Turkey. *J. Radioanal. Nucl. Chem.* **2004**, *259*, 87–92. [CrossRef]

19. Natali, M.; Zanella, A.; Rankovic, A.; Banas, D.; Cantaluppi, C.; Abbadie, L.; Lata, J.C. Assessment of trace metal air pollution in Paris using slurry-TXRF analysis on cemetery mosses. *Environ. Sci. Pollut. Res.* **2016**, *23*, 23496–23510. [CrossRef] [PubMed]

20. Gallego-Cartagena, E.; Morillas, H.; Carrero, J.A.; Madariaga, J.M.; Maguregui, M. Naturally growing grimmiaceae family mosses as passive biomonitors of heavy metals pollution in urban-industrial atmospheres from the Bilbao Metropolitan area. *Chemosphere* **2021**, *263*, 128190. [CrossRef]

21. Cuniasse, B.; Glass, T. Emissions de Gaz à Effet de Serres et de Polluants Atmosphériques en France-Métaux Lourds. CITEPA, Rapport Secten. 2020. https://www.citepa.org/wp-content/uploads/1.3-ML_2020.pdf (accessed on 13 April 2021).

22. IARC. *International Agency for Research on Cancer Monographs on the Evaluation of Carcinogenic Risks to Humans*; IARC: Lyon, France, 2016.

23. Vivanco, M.G.; Gonzalez, M.A.; Palomino, I.; Garrido, J.L.; Querol, X.; Bessagnet, B.; de la Rosa, J.; de la Campa, A.M.S. Modelling arsenic, lead, cadmium and nickel ambient air concentrations in Spain. In Proceedings of the 2011 International Conference on Computational Science and Its Applications, Santander, Spain, 20–23 June 2011; pp. 243–246.

24. Mailler, S.; Menut, L.; Khvorostyanov, D.; Valari, M.; Couvidat, F.; Siour, G.; Turquety, S.; Briant, R.; Tuccella, P.; Bessagnet, B.; et al. CHIMERE-2017: from urban to hemispheric chemistrytransport modeling. *Geosci. Model Dev.* **2017**, *10*, 2397–2423. [CrossRef]

25. Steinnes, E.; Ruhling, A.; Lippo, H.; Makinen, A. Reference materials for large-scale metal deposition surveys. *Accredit. Qual. Assur.* **1997**, *2*, 243–249. [CrossRef]

26. Menut, L.; Bessagnet, B.; Khvorostyanov, D.; Beekmann, M.; Blond, N.; Colette, A.; Coll, I.; Curci, G.; Foret, G.; Hodzic, A.; et al. CHIMERE 2013: A model for regional atmospheric composition modelling. *Geosci. Model Dev.* **2013**, *6*, 981–1028. [CrossRef]

27. Couvidat, F.; Bessagnet, B.; Garcia-Vivanco, M.; Real, E.; Menut, L.; Colette, A. Development of an inorganic and organic aerosol model (Chimere2017β v1.0): Seasonal and spatial evaluation over Europe. *Geosci. Model Dev.* **2018**, *11*, 165–194. [CrossRef]

28. Rouil, L.; Honore, C.; Vautard, R.; Beekmann, M.; Bessagnet, B.; Malherbe, L.; Meleux, F.; Dufour, A.; Elichegaray, C.; Flaud, J.M.; et al. PREV'AIR: An operational forecasting and mapping system for air quality in Europe. *Bull. Am. Meteorol. Soc.* **2009**, *90*, 73–84. [CrossRef]

29. Gareth, J.; Witten, D.; Hastie, T.; Robert, T. *An Introduction to Statistical Learning*; Springer: New York, NY, USA, 2013.

30. Pinheiro, J.C.; Bates, D.M. *Mixed-Effects Models in S and S-Plus*; Springer: Berlin/Heidelberg, Germany, 2000.
31. Demidenko, E. *Mixed Models: Theory and Applications with R*; John Wiley & Sons: Hoboken, NJ, USA, 2013.

32. Cressie, N.; Wickle, C.K. *Statistics for Spatio-Temporal Data*; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2011.

33. Juda-Rezler, K.; Zajusz-Zubek, E.; Reizer, M.; Maciejewska, K.; Kurek, E.; Bulska, E.; Klejnowski, K. Bioavailability of elements in atmospheric PM2.5 during winter episodes at Central Eastern European urban background site. *Atmos. Environ.* **2021**, *245*, 1352–2310. [CrossRef]

34. Guo, G.; Zhang, D.; Wang, Y. Characteristics of heavy metals in size-fractionated atmospheric particulate matters and associated health risk assessment based on the respiratory deposition. *Environ. Geochem. Health* **2021**, *43*, 285–299. [CrossRef]

35. Berg, T.; Røyset, O.; Steinnes, E. Moss (Hylocomium splendens) used as biomonitor of atmospheric trace element deposition: Estimation of uptake efficiencies. *Atmos. Environ.* **1995**, *29*, 353–360. [CrossRef]

36. Berg, T.; Steinnes, E. Use of mosses (Hylocomium splendens and Pleurozium schreberi) as biomonitors of heavy metal deposition: From relative to absolute deposition values. *Environ. Pollut.* **1997**, *98*, 61–71. [CrossRef]

37. Nickel, S.; Schröder, W.; Schaap, M. Estimating atmospheric deposition of heavy metals in Germany using LOTOS-EUROS model calculations and data from biomonitoring programmes. *Pollut. AtmosphÉrique* **2015**, *226*, 2268–3798. [CrossRef]

38. Nickel, S.; Schröder, W.; Wosniok, W.; Harmens, H.; Frontasyeva, M.V.; Alber, R.; Aleksiayenak, J.; Barandovski, L.; Blum, O.; Danielsson, H.; et al. Modelling and mapping heavy metal and nitrogen concentrations in moss in 2010 throughout Europe by applying Random Forests models. *Atmos. Environ.* **2017**, *156*, 146–159. [CrossRef]

39. Shaw, A. *Heavy Metal Tolerance in Plants: Evolutionary Aspects*; CRC Press: Boca Raton, FL, USA, 1990.