# Spatial extrapolation of cadmium concentration in terrestrial mosses using multiple linear regression model predictions across French biogeographical regions

Jérémy Lamouroux, Caroline Meyer, Sébastien Leblond, Isabelle Albert

► To cite this version:

Jérémy Lamouroux, Caroline Meyer, Sébastien Leblond, Isabelle Albert. Spatial extrapolation of cadmium concentration in terrestrial mosses using multiple linear regression model predictions across French biogeographical regions. Environmental Science and Pollution Research, 2025, 10.1007/s11356-025-35985-5. hal-04932902

HAL Id: hal-04932902

https://hal.science/hal-04932902v1

Submitted on 6 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Spatial extrapolation of cadmium concentration in terrestrial mosses using multiple linear regression model predictions across French biogeographical regions

Jérémy Lamouroux<fnref n="1" /> Ⓒ · Caroline Meyer<fnref n="2" /> · Sébastien Leblond<fnref n="2" /> · Isabelle Albert<fnref n="1" />

Received: 21 March 2024 / Accepted: 18 January 2025

© The Author(s) 2025

## Abstract

The French Moss Survey employs forest mosses as indicators to monitor the deposition of atmospheric substances, notably focusing on cadmium (Cd), a known carcinogenic and contributor to respiratory illnesses. This comprehensive study encompasses 55 variables to understand Cd accumulation in terrestrial mosses in France. These variables include moss species, tree cover, biogeographical markers, land use area, proximity to road and rail networks, soil concentration of Cd and atmospheric concentration and deposition of Cd using a physical model. The response variable undergoes a complementary log-log transformation to constrain prediction values within the maximum Cd content in mosses. We have built a regression model to improve predictions, considering the impacts of covariates in France. This model retains biogeographical effects, leading to data segmentation into four distinct biogeographical zones: Atlantic, Continental, Mediterranean and Alpine. Subsequently, zone-specific regression models are explored to refine predictions and consider the impacts of covariates specific to each region, such as those related to railways and roads of the Mediterranean zone. Our biogeographical models effectively mitigate spatial correlation issues and yield accurate predictions, as evidenced by the leave-one-out cross-validation assessment. Compared to ordinary kriging map, the regression prediction maps highlight the contributions of certain covariates, such as the EMEP atmospheric transport model, to areas with high Cd concentrations. Furthermore, these maps exhibit new areas with high (resp. low) Cd concentrations due to high (resp. low) values of the covariates.

**Keywords** Forest mosses · Heavy metal · Linear model · Spatial extrapolation · Prediction map · Biomonitoring

## Introduction

Many metallic elements in ionic or particulate form are present in the atmosphere. These elements become toxic above a certain threshold (Nordberg et al. 2007) and cause major health problems for ecosystems and humans (Gurjar et al. 2010). Among these elements, cadmium (Cd), which is carcinogenic and triggers respiratory disease (Andujar et al. 2010), accumulates in organisms (Zhang and Reynolds 2019). Moreover, Cd is biomagnified in food webs (Butt

et al. 2018) and leads to a loss of biodiversity (Zwart et al. 2010). Most Cd impact studies focus on high concentrations. However, chronic exposure to low concentrations can also increase the risk of cancer (Lequy et al. 2023) and death from heart and respiratory diseases (Lequy et al. 2019). Therefore, assessing the spatial distribution of Cd levels is necessary, especially in rural areas where contamination is considered low compared to urban areas (Vieille et al. 2021).

Few monitoring stations in France measure atmospheric and deposition levels of cadmium (Cd), and those that do are primarily located in urban or peri-urban areas. These measurement sites remain limited due to technical challenges and operating costs. Although these stations meet regulatory requirements for monitoring atmospheric metals, more is needed to provide comprehensive predictions, particularly at a national level where monitoring is sparse. Modelling plays a crucial role in compensating for this lack of coverage. In Europe, for example, cadmium (Cd) deposition and atmospheric concentrations are modelled

Responsible Editor: Elena Maestri

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

<footnote>Jérémy Lamouroux
jeremy.lamouroux@inrae.fr</footnote>

<footnote n="1">¹ MIA-Paris-Mathématiques Et Informatique Appliquées, INRAE-AgroParisTech-Université Paris-Saclay, 22 Place de L'Agronomie, 91477 Palaiseau, France</footnote>

<footnote n="2">² PatriNat (OFB-MNHN), 12 Rue Buffon, 75005 Paris, France</footnote>
by the European Monitoring and Evaluation Programme (EMEP) using an atmospheric transport model (Ilyin et al. 2013; Travnikov and Ilyin 2005), which relies on data about emissions and meteorological conditions from the concerned countries.

Since the 1960s, terrestrial mosses have been used to monitor atmospheric metals, including Cd, at fine scales. This allows for higher sampling densities than conventional deposition measurements (Gusev et al. 2009). Terrestrial mosses morphologies, the absence of roots and a developed vascular system make them excellent sensors of atmospheric contaminants as they depend on atmospheric deposition (Bates 1992). The concentrations of atmospheric Cd are commonly considered proxies for atmospheric deposition. However, it is essential to understand the sources and environmental factors modulating Cd accumulation by mosses (Lequy et al. 2023).

Several studies have demonstrated that Cd concentrations in mosses are influenced by various factors, including Cd deposition modelled by EMEP (Schröder et al. 2019), the percentage of land use determined by the CORINE Land Cover (CLC) (Schröder et al. 2014) and the presence of forests or urban areas (Lequy et al. 2022) or the Cd concentration of the soil on which the mosses grow (Taeprayoon et al. 2023). Moss species and forest cover also play a role in determining metal concentrations (Holy et al. 2009). More importantly, these relationships are contingent. There are three principal approaches to addressing this issue. The first relies on explicative models identifying the most significant covariates explaining metal concentrations in mosses (Lequy et al. 2022). The second approach is based on analyses to understand the links between each component using Spearman correlations, Bartlett's tests and Dunn's tests (Harmens et al. 2012; Boquete et al. 2020). The last approach uses mapping based on regression models or decision trees (Harmens et al. 2011; Schröder et al. 2017). The importance of certain variables becomes apparent in some research works (Lequy et al. 2022), which focused on urban land in France, identifying forest and urban land use as significant variables. Another study on forestry land in France highlighted moss species, atmospheric concentration and deposition as significant variables (Lequy et al. 2017a, b). Interestingly, studies in Germany (Schröder et al. 2019) revealed different significant variables including population density and distance to federal highways. The concept of regionalisation (Vieille et al. 2021) was demonstrated by the varying effects of variables when comparing moss sampling in urban land in two French cities: Paris and Lyon. This regionalisation phenomenon underscored the link between the study area and the relevant covariates. Overall, these assessments contribute to our understanding of the environmental impact of moss accumulation and the distribution of metal concentrations within a region.

This study aims to understand and establish the spatial representation of Cd background levels across mainland France. Cd levels in mosses are influenced by atmospheric deposition, land use and environmental factors such as forest cover, population density and proximity to urban areas or highways. These factors likely vary regionally, reflecting a phenomenon of regionalisation, where different environmental and anthropogenic covariates are significant in distinct areas. Using multiple regression models, we can identify these significant covariates and predict the spatial distribution of Cd in mosses with more precision than traditional methods like ordinary kriging. We first seek to analyse and compare multiple regression models to identify the most relevant environmental and anthropogenic variables affecting Cd concentrations in mosses across France and its biogeographical zones. We then aim to predict the spatial distribution of Cd from these models. To do so, data collected from 445 sites of the BRAMM (Biosurveillance des Retombées Atmosphériques Métalliques par les Mousses) network were analysed.

# Materials and methods

## Measurements Cd in forest mosses across France

### Study sites and moss sampling

During the BRAMM 2016 campaign, moss samples were collected from 445 forested sites evenly spread across France. These sites represent background pollution levels. Moss samples were collected following the standardised protocol of ICP Vegetation (Monitoring ICP Manual, 2020). In brief, a single moss species was collected from areas distant from any known local pollution sources. Each site covered an area of 2500 m² and was located below tree canopies. One of five different moss species was collected for each site to obtain samples from 445 sites: *Hypnum cupressiforme* Hedw., 1801 (Hc, 227 samples), *Pseudoscleropodium purum* (Hedw.) M.Fleisch., 1923 (Pp, 186 samples), *Thuidium tamariscinum* (Hedw.) Schimp., 1852 (Tt, 27 samples), *Pleurozium schreberi* (Willd. ex Brid.) Mitt., 1869 (Ps, 4 samples), and *Hylocomium splendens* (Hedw.) Schimp., 1852 (Hs, 1 sample) (Fig. 1a), respectively named Hc, Pp, Tt, Ps and Hs hereafter. The five species could be distinguished in the field using morphological criteria that could be observed with a magnifying glass (magnification ×10; lens diameter 21 mm). For each sample, 10 fragments of the same moss species carpets were collected. *Hypnum cupressiforme* was always collected on dead wood, while other species were collected on humus (composed solely of decomposing organic matter) or soil (composed of a mixture of inorganic and organic matter). Each sample was stored in plastic and brought to the

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 1** a The 445 observed forest sites (BRAMM sites), depending on the moss species. b The four biogeographical zones, defined by the official boundaries established in the REGBIOFR layer (Roekaerts 2012). Hc, *Hypnum cupressiforme* Hedw.; Pp, *Pseudoscleropodium*</figcaption></figure>

purum (Hedw.) M.Fleisch.; Ps, *Pleurozium schreberi* (Willd. ex Brid.) Mitt.; Tt, *Thuidium tamariscinum* (Hedw.) Schimp.; Hs, *Hylocomium splendens* (Hedw.) Schimp

laboratory within the same week. At the laboratory, samples were carefully cleaned from all dead material and attached litter. The whole moss shoot was kept. Samples were dried at room temperature before grinding with an automatic non-polluting titanium grinder (Pulverisette®14, Fritsch, Germany). All necessary precautions were taken to avoid contamination while working with trace levels, including using powder-free gloves, plastic and ceramic pincers and scissors.

## Analysis of the collected moss samples and quality assurance

A total of 250 mg of homogeneous powder was dried for 48 h at 24 °C and sent to USRAVE (INRAE-Centre de Bordeaux) for analysis. The powder was mineralised with a mixture of HF/HNO₃/H₂O₂ (concentrations expressed at 103 °C after accounting for water loss on a subsample). Cadmium in mosses was analysed by inductively coupled plasma mass spectrometry (ICP-MS) by the USRAVE (INRAE-Centre de Bordeaux) (limit of quantification was 0.02 µg/g; analytical uncertainty is 20%). The device used by the USRAVE was an Agilent 7700 × spectrometer.

The ICP-MS was calibrated by external calibration, using multi-element solutions produced every 6 months. The internal calibration procedure involved the addition of the standard. USRAVE uses the NIST 1547 (peach leaves) to validate their results as an internal Certified Reference Material for trace elements (CRM). We had 3 other CRMs analysed namely, NIST 1573a (tomato leaves), as well as European moss survey programme reference materials—Moss M2 and Moss M3—were processed with each batch of samples. The

Moss M2 and Moss M3 standards were powders of *Pleurozium schreberi* from two Finnish sites (Berg et al. 1997).

## Variables of interest

Table 1 presents the 55 constructed variables from information collected at each site, which could explain the variability of Cd concentration in mosses. In this study, the substrate from which the moss was taken was not retained as a variable due to a bias in the protocol as the sampling method created a strong link between substrate and moss species. These 55 variables correspond to 65 covariates, which can be categorised as follows:

i). Categorical covariates encompass 3 variables which are the following: the moss species collected corresponding to 5 covariates (Hc, Hs, Pp, Ps and Tt), the biogeographical zones corresponding to 4 covariates (Atlantic, continental, Mediterranean and Alpine), and the different tree canopies on the sampling site corresponding to 4 covariates (deciduous, coniferous, mixed deciduous tree and mixed coniferous tree).

ii). A quantitative covariate encompasses 1 variable which corresponds to 1 covariate, the altitude of the sampling location (in metres).

iii). Quantitative covariates involve 36 variables, i.e. 36 covariates, land use percentage covariates, expressed as the percentage of land use area within buffer zones of 1 km, 5 km, 10 km and 15 km from the sampling sites, calculated using the 2018 CLC database (with a resolution of 100 × 100 m).
<table><caption>**Table 1** List of 55 constructed variables and their description as potential predictor of Cd accumulation in mosses in 2016</caption><thead><tr><th>Variable</th><th>Code</th><th>Source</th><th>Calculating variables with CLC codes</th><th>Number and name of covariates<sup>a</sup></th></tr></thead><tbody><tr><td>Moss species</td><td>Moss</td><td>BRAMM<sup>a</sup></td><td></td><td>5: Hc, Pp, Ps, Tt, Hs</td></tr><tr><td>Tree cover</td><td>Tree</td><td>BRAMM<sup>a</sup></td><td></td><td>4: deciduous, coniferous, mixed deciduous tree, mixed coniferous tree</td></tr><tr><td>Biogeographical zone</td><td>Biogeo</td><td>INPN<sup>b</sup></td><td></td><td>4: Continental, Atlantic, Mediterranean, Alpine</td></tr><tr><td>Altitude</td><td>Altitude</td><td>Copernicus<sup>c</sup></td><td></td><td>1</td></tr><tr><td>Forested land uses (1, 5, 10, and 15 km radius)</td><td>Forest_*</td><td>CLC<sup>d</sup></td><td>244 + 311 + 312 + 313 + 323 + 324</td><td>4</td></tr><tr><td>Urban land uses (1, 5, 10, and 15 km radius)</td><td>Urban_*</td><td>CLC<sup>d</sup></td><td>111 + 112 + 122 + 141 + 142</td><td>4</td></tr><tr><td>Industrial land uses (1, 5, 10, and 15 km radius)</td><td>Industrial_*</td><td>CLC<sup>d</sup></td><td>121 + 123 + 124 + 131 + 132 + 133</td><td>4</td></tr><tr><td>Pasture land uses (1, 5, 10, and 15 km radius)</td><td>Pasture_*</td><td>CLC<sup>d</sup></td><td>231</td><td>4</td></tr><tr><td>Agricultural land uses (1, 5, 10, and 15 km radius)</td><td>Agricultural_*</td><td>CLC<sup>d</sup></td><td>211 + 212 + 213 + 241 + 242 + 243</td><td>4</td></tr><tr><td>Fruit land uses (1, 5, 10, and 15 km radius)</td><td>Fruits_*</td><td>CLC<sup>d</sup></td><td>221 + 222 + 223</td><td>4</td></tr><tr><td>Vegetation land uses (1, 5, 10, and 15 km radius)</td><td>Vegetation_*</td><td>CLC<sup>d</sup></td><td>321 + 322 + 333</td><td>4</td></tr><tr><td>Water land uses (1, 5, 10, and 15 km radius)</td><td>Water_*</td><td>CLC<sup>d</sup></td><td>335 + 411 + 412 + 511 + 512</td><td>4</td></tr><tr><td>Sea land uses (1, 5, 10, and 15 km radius)</td><td>Sea_*</td><td>CLC<sup>d</sup></td><td>331 + 421 + 422 + 423 + 521 + 522 + 523</td><td>4</td></tr><tr><td>Railway (1, 5, 10, and 15 km radius)</td><td>Railway_*</td><td>CLC<sup>d</sup></td><td></td><td>4</td></tr><tr><td>Road (1, 5, 10, and 15 km radius)</td><td>Road_*</td><td>CLC<sup>d</sup></td><td></td><td>4</td></tr><tr><td>Cd air deposition (1, 3, 6, 9, and 12 month)</td><td>EMEPDep<sup>+</sup></td><td>EMEP<sup>e</sup></td><td></td><td>5</td></tr><tr><td>Cd air concentration</td><td>EMEP-air</td><td>EMEP<sup>e</sup></td><td></td><td>1</td></tr><tr><td>Cd soil total</td><td>RMQS<sub>tot</sub></td><td>RMQS<sup>f</sup></td><td></td><td>1</td></tr></tbody></table>

iv). Quantitative covariates also involve 8 variables, i.e. 8 covariates, distance-related covariates, which measure the distances in metres to roads and railways within the same buffer zones (1 km, 5 km, 10 km and 15 km).

v). Two quantitative covariates derived from the EMEP model include two separate measures represented by 6 variables. Firstly, the atmospheric deposition of Cd in (g/km²/year) was calculated by summing deposits over the 1, 3, 6, 9 or 12 months preceding the collection

period, moving backwards to derive 5 variables, i.e., 5 covariates (EMEP\_dep1, EMEP\_dep3, EMEP\_dep6, EMEP\_dep9, EMEP\_dep12). Secondly, the atmospheric concentration of Cd in (µg/m³/year), modelled in 2016 with a resolution of 0.1°×0.1°, corresponds to 1 variable, i.e., 1 covariate.

vi). A quantitative covariate, represented by 1 variable, i.e. 1 covariate, of Cd soil concentration value (mg/kg²)

<footnote><sup>a</sup>Biosurveillance des Retombées Atmosphériques Métalliques par les Mousses, https://bramm.mnhn.fr</footnote>

<footnote><sup>b</sup>L'Inventaire National du Patrimoine Naturel, https://inpn.mnhn.fr</footnote>

<footnote><sup>c</sup>Copernicus, https://www.copernicus.eu/fr</footnote>

<footnote><sup>d</sup>Corine Land Cover, https://land.copernicus.eu/pan-european/corine-land-cover</footnote>

<footnote><sup>e</sup>European Monitoring and Evaluation Programme, https://www.emep.int/</footnote>

<footnote><sup>f</sup>Réseau de Mesures de la Qualité des Sols, https://www.gissol.fr/le-gis/programmes/rmqs-34</footnote>

<footnote>*The length of radius in km</footnote>

<footnote>+The number of months summed</footnote>

<footnote>The total of covariates 65</footnote>
modelled in 2018 for the period 2000–2010 by the RMQS network (Gis Sol) in a grid cell of 2×2 km.

## Exploration of biogeographical zones

France can be delineated into various zones based on regional, departmental, city or economic criteria. The primary consideration in such divisions is often the geographical extent. In our study, we adopted the division of Metropolitan France into four biogeographical zones, as defined in the Habitats Directive (92/43/EEC). We used the REGBI-OFR layer which defined the official boundaries of biogeographical regions (Fig. 1b). These divisions, endorsed by the Habitats Committee in 2012 (Roekaerts 2012), provided a systematic framework for understanding and managing biodiversity within each region.

The division of France into biogeographical zones was a crucial aspect of our study methodology. This regionalized approach enhanced our ability to capture nuanced environmental factors and their influence on metal accumulation within mosses across Metropolitan France.

## Multiple regression model

### Complementary log-log link

Multiple regression models are widely utilised in data analysis, assuming a linear relationship between response variables, such as Cd concentration in mosses, and explanatory covariates, as described in the “Variables of interest” section. These models offer interpretability through coefficients, indicating the strength and direction of relationships. Linear models also help identify outliers and influential observations. To achieve a more normal distribution of the response variable and align predictions with realistic thresholds, we explored transforming the response variable.

Our method used a complementary log-log (CLL) linear model; the inverse of the Gumbel cumulative distribution function yields the CLL link function (Shim et al. 2023). The model is the following:

$$ln(ln(\frac{Cd_{max}}{Y_i})) = \alpha + \sum_{j=1}^{N} \beta_j X_i^j + \epsilon_i,$$

where ln(.) is the natural logarithm function, $Y_i$ is the Cd concentration in mosses i in (µg/g), $Cd_{max}$ fixed to 1.5 µg/g according to the maximal expected value of Cd in mosses, $\alpha$ and $\beta_j$ are the intercept and slope regression coefficients, respectively, $\sum_{j=1}^{N} X_i^j$ are the sum of N independent covariates $X_i^j$ of interest and $\epsilon_i \sim N(0, \sigma^2)$, where $N(0, \sigma^2)$ is a normal distribution with a mean of 0 and a variability of $\sigma^2$.

## Covariate selection

In the initial step, we filtered out 55 variables (excluding the 3 qualitative variables) with more than 95% of zero values within each biogeographical zone. For example, 95% of the 445 observations had a 0 value for the covariate Marin in the Continental zone. This ensured that the indication provided by the remaining covariates was rare and pertinent, as an overwhelming presence of zeros might render the information irrelevant. Consequently, the number of covariates decreased for each zone. For instance, in the case of the regression model encompassing France, the covariate ‘Fruits’ with a buffer of a 1-km radius was removed due to this criterion.

Following this initial step, we employed the Spearman correlation to investigate collinearity problems between covariates. This non-parametric measure assesses the strength and direction of the relationship between two ranked variables. Spearman correlation values range from -1 to +1, with +1 indicating a perfect positive correlation (Schober et al. 2018). If the Spearman correlation between multiple variables exceeded 0.75, we addressed this by considering two options. Firstly, if the same covariates with different buffers exhibit high correlation (e.g. ‘Forest_1’, ‘Forest_5’, ‘Forest_10’ and ‘Forest_15’), we retained the one with the smallest buffer size (in this case, ‘Forest_1’). Secondly, if there was a correlation between two distinct covariates, such as ‘Forest_1’ and ‘Culture_1’, and the correlation value exceeded 0.75, we retained the covariate with a lower correlation and a separate buffer, here ‘Forest_1’.

To produce the actual selection of variables, i.e. to retain in the model only those variables that influence the response, the Cd concentration in the mosses, we employed a backward-forward selection approach using the Akaike information criterion (AIC), balancing goodness of fit and model complexity. Minimising the AIC value ensured the selection of the most appropriate regression model, effectively reducing the dimensionality while retaining the most influential explicative covariates for our response variable.

## Residual spatial correlation

After selecting the best model, we aimed at investigating whether the spatial correlation of observations remained when covariates are not considered. Unaccounted environmental factors may have contributed to similarities between observations. Thus, determining whether this spatial correlation persisted after introducing selected covariates was essential. The Spatial Correlation Moran Index (Banerjee et al. 2003) is a metric designed to assess whether similar values tend to cluster or disperse across a geographic area and can be used for this purpose. Positive values of the Moran Index, calculated on the response or residuals of the multiple
regression, indicate clustering. This suggests that similar values are grouped together. In contrast, negative values imply dispersion, indicating that similar values are spread apart. A value of 0 suggests the absence of a spatial pattern.

## Prediction maps

### Drawing up concentration maps

It is important to extrapolate regression model predictions to comprehensively assess continuous and fine-tuned territorial variations in Cd concentrations across France. We established a spatial grid covering France with a resolution of 2×2 km to match the size of the moss sampling area, which was a circular area of 2500 m². Predictors were calculated by inverting the multiple regression CLL linear models, providing predictions ŷ(s) (expressed as Hc mosses under deciduous cover tree) at each point (s) of the spatial grid, totalling 139,673 points across France. This calculation required obtaining the values of each covariate retained in the models, sourced from relevant databases covering the territory.

We also produced a map of the standard deviations of each prediction (obtained using the predict function from terra R package, which approximates this variance using the delta method). This was useful for visualising the uncertainty associated with the calculated prediction.

Also, to assess the contribution of our covariates to our predictions, we compared our prediction maps (point prediction and standard deviation) with maps obtained using the classic ordinary kriging (OK) method, which relied solely on observations (response) and their spatial covariance to obtain predictions (Cressie 1990). The observations underwent a logarithmic transformation to improve normality.

## Investigation of the predictive performance

The precision of the extrapolation methods was assessed using the leave-one-out cross-validation (LOOCV) technique (Wilcox et al. 2012). LOOCV is a widely used approach for comparing the accuracy of extrapolation methods. It involves systematically excluding each sampling site. Subsequently, the value at the excluded site was estimated based on the prediction at the site calculated with the model excluding the site. This iterative process was repeated until all sampling sites have prediction values, resembling an extreme case of k-fold cross-validation, where k equals the total data points in the dataset. Test points were used to calculate the root mean squared error (RMSE), which quantified the average difference between the predictions of the dataset and the observed values. This measurement helped assess extrapolation performance (Stone 1974; Tyagi et al. 2022):

$$LOOCV RMSE = \sqrt{\sum_{i=1}^{M} \frac{(y_i - \hat{y}_i)^2}{M}},$$

where at site i, $y_i$ is the Cd concentration observed, $\hat{y}_i$ is the concentration prediction fitted without $y_i$ and M is the total number of observations in the model. The lower the LOOCV value, the better the model predictions match the observed data.

All analyses were conducted using R version 4.2.2. Spatial data manipulation and analysis were performed using the ‘terra’ package (version 1.7–39). Linear regression models were fitted using the ‘lm’ function, and model selection was based on the Akaike information criterion (AIC), calculated using the ‘AIC’ function. Predictions were generated using the ‘predict’ function, and spatial visualisations were created using the ‘ggplot’ function.

## Results

### Cadmium concentration in CRMs and BRAMM mosses

The assessment of our 3 CRMs (Table 2) from the 15 series of analyses carried out during the 2016 campaign validated the accuracy of our results. There was no difference between the mean Cd concentration, including its standard deviation, associated with the reference values and those obtained from our analyses.

### Model and map encompassing France

#### Model selection

We started by analysing a model encompassing France, utilising a dataset of 445 observations and 55 variables corresponding to 65 covariates. We removed variables with excessive zeros and high collinearity, as detailed in the “Materials and methods” section. Lastly, we retained 27 variables (Table 3) corresponding to 37 covariates.

The AIC backward-forward selection retains 9 relevant variables for modelling France (Table 3): moss species, tree cover, biogeographical zone, urban land uses with a buffer of 5 km, industrial land uses with a buffer of 5 km, agricultural land uses with a buffer of 1 km, sea land uses with a buffer of 10 km, Cd soil total RMQS and Cd air concentration EMEP. The selected model reaches an AIC value of 167.81. Only 7 of these covariates demonstrate a significant p-value < 0.05, as shown in Table 4.

Three significant qualitative covariates influence the concentration of Cd in mosses: moss species, tree cover and biogeographical zone. These covariates highlight the
<table><caption>**Table 2** Elementary statistics of Cd concentrations analysed in reference materials (SRM 1573a, Moss M2, and Moss M3, n=15) and in 445 moss samples taken during the 2016 BRAMM campaign</caption><thead><tr><th rowspan="2"></th><th colspan="3">References materials **</th><th rowspan="2">Cd<br>(BRAMM<br>2016)</th></tr><tr><th>Moss M2</th><th>Moss M3</th><th>NIST 1573a</th></tr></thead><tbody><tr><td>Number of values</td><td>15</td><td>15</td><td>15</td><td>445</td></tr><tr><td>Minimum</td><td>0.433</td><td>0.0931</td><td>1.350</td><td>0.029</td></tr><tr><td>1st quartile</td><td>0.453</td><td>0.107</td><td>1.440</td><td>0.099</td></tr><tr><td>Median</td><td>0.460</td><td>0.111</td><td>1.470</td><td>0.135</td></tr><tr><td>3rd quartile</td><td>0.477</td><td>0.113</td><td>1.505</td><td>0.199</td></tr><tr><td>Maximum</td><td>0.558</td><td>0.127</td><td>1.640</td><td>1.160</td></tr><tr><td>Mean</td><td>0.468</td><td>0.111</td><td>1.473</td><td>0.174</td></tr><tr><td>Standard deviation</td><td>0.028</td><td>0.009</td><td>0.079</td><td>0.138</td></tr><tr><td>Coefficient of variation (%)</td><td>6</td><td>8</td><td>5</td><td>79</td></tr><tr><td>Ratio max/min</td><td>1</td><td>1</td><td>1</td><td>40</td></tr><tr><td>Atypical high threshold*</td><td>0.513</td><td>0.123</td><td>1.602</td><td>0.349</td></tr><tr><td>Number of typical threshold values</td><td>1</td><td>2</td><td>1</td><td>28</td></tr></tbody></table>

importance of moss species, the specific forest cover where the mosses grow and the biogeographical zone in which they are situated. Two quantitative covariates exhibit high significance levels (*p* < 0.001): Sea_10, with a negative effect, and EMEP-air, with a positive effect. Conversely, EMEP-air indicates that as Cd air concentration increases, Cd concentrations in mosses also rise. On the contrary, the percentage of sea land uses around the moss collection area are associated with lower Cd concentrations in mosses. Another quantitative covariate, Agricultural_1, shows statistical significance at *p* < 0.01, with a negative effect. Higher percentages of agricultural land uses around the moss collection area are associated with lower Cd concentrations in mosses. Finally, the last significant covariate, RMQS_tot_1, reaches statistical significance at *p* < 0.05, and its influence is positively oriented (Table 4). When the concentration of Cd in the soil is high at the location where the moss has been collected, the concentration of Cd in the moss is higher. Urban_5 and Industrial_5 are non-significant but remain in the model because they reduce the AIC and improve the model's goodness of fit compared to a model without them.

The application of Moran's index, as described in the "Material and methods" section, assesses whether the model with 9 variables fully accounts for spatial effects or if residual spatial autocorrelation remains after accounting for these covariates. Table 5 indicates that residual spatial autocorrelation persists, suggesting that the covariates do not fully explain the spatial patterns observed. Note that Moran's index value is divided by 2 in our analysis.

## Prediction maps

We constructed a prediction map of Cd concentration throughout France using the 9 retained variables. Figure 2a and b depict the extrapolation of Cd concentration with moss Hc under deciduous forest cover, along with standard deviation maps. The white borders delineating biogeographical zones indicate the relevance of the variable to the model. The standard deviation maps illustrate prediction accuracy, with higher or lower predictions corresponding to higher standard deviations. Supplementary data provide additional maps for other moss species and forest cover types. Figure 3a displays the ordinary kriging of Cd concentration in mosses, while Fig. 3b represents the standard deviation of ordinary kriging values.

Our spatial analysis reveals elevated predictions of Cd concentrations along the border with Belgium and in northern France, particularly in proximity to Paris. Another zone with higher Cd concentration predictions has also been identified near Marseille in southern France. Conversely, within the Mediterranean zone, a distinct low prediction of Cd concentrations is observed. The standard deviation map highlights significant variability, with pronounced high- and low-concentration prediction fluctuations. Notably, the southern region of France exhibits substantial variability, corresponding to areas of elevated predictions. A similar pattern is observed near Paris.

It is interesting to compare the previous map of the model's prediction across France with the ordinary kriging of Cd concentration in mosses (Fig. 3a). Both indicate

<footnote>*Atypical high threshold = 3rd quart. + 1,5 * (3rd quart. - 1st quart.)</footnote>

<footnote>**Values of reference materials: Moss M2=0.454±0.019 µg/g; Moss M3=0.106±0.005 µg/g; NIST 1573a=1.52±0.04 µg/g</footnote>
<figure><img src="image_5.png" /><figcaption>**Table 3** List of variables before and after the backwards-forwards Akaike index criterion (AIC) selection on the complementary log-log linear model in France and each biogeographical zone. Full model = variables selected after removing those for which more than 95% of the values were equal to 0, then after selection by Spearman correlation; relevant = variables selected after backward-forward AIC selection; total of variables = total number of variables selected with the number of significant variables in brackets (p < 0.05)</figcaption></figure>

<table><thead><tr><th rowspan="3">Variable</th><th rowspan="3">Code</th><th colspan="10">Selection</th></tr><tr><th colspan="2">France (n=445)</th><th colspan="2">Atlantic zone (n=220)</th><th colspan="2">Continental zone (n=147)</th><th colspan="2">Mediterranean zone (n=49)</th><th colspan="2">Alpine zone (n=29)</th></tr><tr><th>Full model</th><th>Relevant</th><th>Full model</th><th>Relevant</th><th>Full model</th><th>Relevant</th><th>Full model</th><th>Relevant</th><th>Full model</th><th>Relevant</th></tr></thead><tbody><tr><td>Moss species</td><td>Moss<sup>a</sup></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Tree cover</td><td>Tree<sup>a</sup></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Biogeographical zone</td><td>Bioge<sup>b</sup></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Altitude</td><td>Altitude<sup>c</sup></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Forested land uses (1, 5, 10, and 15 km radius)</td><td>Forest_<sup>nd</sup></td><td>1,5</td><td></td><td>1,5</td><td>1</td><td>1,5</td><td>1,5</td><td>1,5</td><td>1,5</td><td>1,5</td><td>1,5</td></tr><tr><td>Urban land uses (1, 5, 10, and 15 km radius)</td><td>Urban_<sup>nd</sup></td><td>1,5</td><td>5</td><td>1,5</td><td>5</td><td>1,5</td><td>5</td><td>1,5</td><td>1,5</td><td>1,5</td><td>1,5</td></tr><tr><td>Industrial land uses (1, 5, 10, and 15 km radius)</td><td>Industrial_<sup>nd</sup></td><td>5,10</td><td>5</td><td>5,10</td><td>5</td><td>5</td><td>5</td><td>5</td><td>5,10</td><td>5,10</td><td>5,10</td></tr><tr><td>Pasture land uses (1, 5, 10, and 15 km radius)</td><td>Pasture_<sup>nd</sup></td><td>1,5</td><td></td><td>1,5</td><td>1,5</td><td>1,5</td><td>5</td><td>1,5</td><td>1</td><td>1,5</td><td>1,5</td></tr><tr><td>Agricultural land uses (1, 5, 10, and 15 km radius)</td><td>Agricultural_<sup>nd</sup></td><td>1,5</td><td>1</td><td>5</td><td>1,5</td><td>1,5</td><td>5</td><td>1,5</td><td>1,5</td><td>1,5</td><td>1,5</td></tr><tr><td>Fruits land uses (1, 5, 10, and 15 km radius)</td><td>Fruits_<sup>nd</sup></td><td>5,10</td><td></td><td>5,10</td><td>5,10</td><td>1,5</td><td>5</td><td>1,5</td><td>5</td><td>10</td><td>10</td></tr><tr><td>Vegetation land uses (1, 5, 10, and 15 km radius)</td><td>Vegetation_<sup>nd</sup></td><td>1,5</td><td></td><td>5,10</td><td>1,5</td><td>5</td><td>1,5</td><td>1,5</td><td>1,5</td><td>1,5</td><td>1,5</td></tr><tr><td>Water land uses (1, 5, 10, and 15 km radius)</td><td>Water_<sup>nd</sup></td><td>5,10</td><td></td><td>5,10</td><td>5,10</td><td>5,10</td><td>5</td><td>5,10</td><td>5,10</td><td>5,10</td><td>5,10</td></tr><tr><td>Sea land uses (1, 5, 10, and 15 km radius)</td><td>Sea_<sup>nd</sup></td><td>10</td><td>10</td><td>5,10</td><td>10</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Railway (1, 5, 10, and 15 km radius)</td><td>Railway_<sup>nd</sup></td><td>5</td><td></td><td>5,10</td><td></td><td>5</td><td>5</td><td>5,10</td><td>5,10</td><td>5</td><td>5</td></tr><tr><td>Road (1, 5, 10, and 15 km radius)</td><td>Road_<sup>nd</sup></td><td>5</td><td></td><td>5</td><td></td><td>5</td><td>5</td><td>5,10</td><td>5,10</td><td>5,10</td><td>5,10</td></tr><tr><td>Cd air deposition (1, 3, 6, 9, and 12 month)</td><td>EMEP_<sup>dep+e</sup></td><td>1,3</td><td></td><td>1,3</td><td></td><td>1</td><td>1</td><td>1,3</td><td></td><td>1</td><td>1</td></tr><tr><td>Cd air concentration</td><td>EMEP_<sup>air+e</sup></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Cd soil total</td><td>RMQS_<sup>totf</sup></td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Total of variables</td><td></td><td>27</td><td>9 (7*)</td><td>27</td><td>10 (8*)</td><td>23</td><td>12 (5*)</td><td>29</td><td>20 (11*)</td><td>22</td><td>22 (3*)</td></tr><tr><td>AIC</td><td></td><td>167.81</td><td></td><td>12.79</td><td></td><td>87.90</td><td></td><td>-9.47</td><td></td><td>181.27</td><td></td></tr></tbody></table>

<footnote><sup>a</sup>Biosurveillance des Retombées Atmosphériques Métalliques par les Mousses, https://bramm.mnhn.fr/</footnote>

<footnote><sup>b</sup>L'Inventaire National du Patrimoine Nature, https://inpn.mnhn.fr</footnote>

<footnote><sup>c</sup>Copernicus, https://www.copernicus.eu/fr</footnote>

<footnote><sup>d</sup>Corine Land Cover, https://land.copernicus.eu/pan-european/corine-land-cover</footnote>

<footnote><sup>e</sup>European Monitoring and Evaluation Programme, https://www.emep.int/</footnote>

<footnote><sup>f</sup>Réseau de Mesures de la Qualité des Sols, https://www.gissol.fr/le-gis/programmes/rnqs-34</footnote>

<footnote><sup>g</sup>The length of radius in km</footnote>

<footnote><sup>h</sup>The number of months summed</footnote>

<footnote><sup>i</sup>Covariates with significant p-value < 0.05</footnote>
<caption>**Table 4** Significant covariates ($p < 0.05$) and their effect (negative or positive) on the model in France and each biogeographical zone. Moss Hs, Moss Pp, Moss Ps, and Moss Tt=covariates associated with the variable moss species and expressed relatively to the Hc moss species (Hs, *Hylocomium splendens*; Pp, *Pseudosclopodium purum*; Ps, *Pleurozium schreberi*; Tt, *Thuidium tamariscinum*; Hc, *Hypnum cupressiforme*). Biogeo ATL, Biogeo CON,</caption>

Biogeo MED=covariates associated with the variable Biogeographical zone and expressed relatively to the Alpine zone (ATL, Atlantic zone; CON, Continental zone; MED, Mediterranean zone). Coniferous trees, mixed coniferous tree, mixed deciduous tree=covariates associated with the variable tree cover and expressed relatively to the deciduous trees

<table><thead><tr><th rowspan="2">Code</th><th colspan="2">France</th><th colspan="2">Atlantic zone</th><th colspan="2">Continental zone</th><th colspan="2">Mediterranean zone</th><th colspan="2">Alpine zone</th></tr><tr><th>Effect</th><th>p_value</th><th>Effect</th><th>p_value</th><th>Effect</th><th>p_value</th><th>Effect</th><th>p_value</th><th>Effect</th><th>p_value</th></tr></thead><tbody><tr><td>Intercept</td><td>-</td><td>&lt;0.001</td><td>-</td><td>&lt;0.001</td><td>+</td><td>&lt;0.05</td><td></td><td></td><td></td><td></td></tr><tr><td>Moss Hs</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Moss Pp</td><td>-</td><td>&lt;0.001</td><td>-</td><td>&lt;0.001</td><td>+</td><td>&lt;0.001</td><td></td><td></td><td></td><td></td></tr><tr><td>Moss Ps</td><td>-</td><td>&lt;0.05</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Moss Tt</td><td>-</td><td>&lt;0.001</td><td>-</td><td>&lt;0.01</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Biogeo ATL</td><td>+</td><td>&lt;0.05</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Biogeo CON</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Biogeo MED</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Coniferous trees</td><td>+</td><td>&lt;0.001</td><td>+</td><td>&lt;0.05</td><td>+</td><td>&lt;0.001</td><td></td><td></td><td></td><td></td></tr><tr><td>Mixed coniferous tree</td><td>+</td><td>&lt;0.001</td><td>+</td><td>&lt;0.001</td><td>+</td><td>&lt;0.05</td><td>-</td><td>&lt;0.05</td><td></td><td></td></tr><tr><td>Mixed deciduous tree</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Altitude</td><td></td><td></td><td></td><td></td><td>-</td><td>&lt;0.05</td><td></td><td></td><td></td><td></td></tr><tr><td>Forest_1</td><td></td><td></td><td>+</td><td>&lt;0.05</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Forest_5</td><td></td><td></td><td></td><td></td><td></td><td></td><td>-</td><td>&lt;0.05</td><td></td><td></td></tr><tr><td>Urban_1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Urban_5</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Industrial_5</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Industrial_10</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Pasture_1</td><td></td><td></td><td></td><td></td><td></td><td></td><td>+</td><td>&lt;0.01</td><td></td><td></td></tr><tr><td>Pasture_5</td><td></td><td></td><td>-</td><td>&lt;0.05</td><td>+</td><td>0.05</td><td></td><td></td><td></td><td></td></tr><tr><td>Agricultural_1</td><td>-</td><td>0.01</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Agricultural_5</td><td></td><td></td><td></td><td></td><td></td><td></td><td>-</td><td>0.05</td><td></td><td></td></tr><tr><td>Fruits_5</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Fruits_10</td><td></td><td></td><td>+</td><td>&lt;0.05</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Vegetation_1</td><td></td><td></td><td></td><td></td><td></td><td></td><td>+</td><td>&lt;0.01</td><td></td><td></td></tr><tr><td>Vegetation_5</td><td></td><td></td><td></td><td></td><td>+</td><td>&lt;0.05</td><td>-</td><td>&lt;0.01</td><td></td><td></td></tr><tr><td>Water_10</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Water_15</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Sea_5</td><td></td><td></td><td></td><td></td><td></td><td></td><td>-</td><td>&lt;0.05</td><td></td><td></td></tr><tr><td>Sea_10</td><td>-</td><td>&lt;0.001</td><td>-</td><td>&lt;0.05</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Railway_5</td><td></td><td></td><td></td><td></td><td></td><td></td><td>+</td><td>&lt;0.01</td><td></td><td></td></tr><tr><td>Railway_10</td><td></td><td></td><td></td><td></td><td></td><td></td><td>-</td><td>&lt;0.05</td><td></td><td></td></tr><tr><td>Road_5</td><td></td><td></td><td></td><td></td><td></td><td></td><td>-</td><td>&lt;0.01</td><td></td><td></td></tr><tr><td>Road_10</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>EMEPDep1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>-</td><td>0.05</td></tr><tr><td>EMEP-air</td><td>+</td><td>&lt;0.001</td><td>+</td><td>&lt;0.001</td><td></td><td></td><td>+</td><td>&lt;0.05</td><td>+</td><td>0.05</td></tr><tr><td>RMQS Tot</td><td>-</td><td>&lt;0.05</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

higher concentrations in mosses near Lille and Paris, but an important difference is observed near Marseille. The regression models' predictions are notably higher near Marseille and are related to the highest values modelled by EMEP. The regression model's prediction map also reveals a lower

coastline contamination level attributed to the sea covariate within a 10-km buffer. The maps covering France (Fig. 2a and Fig. 3a) highlight the possibility of high Cd concentrations in mosses near the Belgian border and the city of Lille. However, the standard deviation of these predictions
<table><caption>**Table 5** Moran's index (I-statistics) before and after considering covariates in the models. Atlantic, Continental, Mediterranean, Alpine = biogeographical zone. NS is not significant (p > 0.05)</caption><thead><tr><th rowspan="2"></th><th colspan="2">Without considering covariates</th><th colspan="2">Considering covariates</th></tr><tr><th>I-statistics</th><th>p-value</th><th>I-statistics</th><th>p-value</th></tr></thead><tbody><tr><td>France</td><td>0.307</td><td>p < 0.001</td><td>0.154</td><td>p < 0.05</td></tr><tr><td>Atlantic zone</td><td>0.424</td><td>p < 0.001</td><td>0.011</td><td>p < 0.05</td></tr><tr><td>Continental zone</td><td>0.091</td><td>p < 0.001</td><td>-0.010</td><td>NS</td></tr><tr><td>Mediterranean zone</td><td>0.064</td><td>p < 0.001</td><td>0.201</td><td>NS</td></tr><tr><td>Alpine zone</td><td>0.476</td><td>p < 0.001</td><td>-0.181</td><td>NS</td></tr></tbody></table>

is high (Fig. 2b and Fig. 3b). The standard deviation map shows greater variability where the high and low concentration predictions lie.

Table 6 gives the LOOCV of the 9-variable model over France (0.296) calculated zone by zone. Additionally, the LOOCV for ordinary kriging is reported as 0.491 over France. Comparing the two methods regarding the LOOCV criterion (Table 6), our maps demonstrate predictions closer to the observations, instilling confidence in the other's predictions in areas lacking observation data for comparison.

## Models and maps by biogeographical zone

Given the statistical significance of the biogeographical zone covariate and the persistence of spatial residuals, we explore more refined models within each biogeographical zone.

## Model selection

Table 3 presents the full model for each biogeographical zone, including the number of observations and covariates in the model after removing variables with too many zeros and those that are too collinear. Additionally, Table 3 displays the relevant model selected based on the smallest AIC criterion and the significant covariates in these models. Table 4 provides the significance levels of the covariates.

In the Atlantic zone, a model built with 220 observations using backward-forward AIC selection retains 10 relevant variables, including 8 with significant associations with the response variable. Two qualitative variables, moss species and tree cover, are retained due to significant differences from the intercept (Hc mosses and deciduous trees). Pp mosses (p < 0.001) and Tt mosses (p < 0.01) exhibit significant differences from Hc mosses under deciduous trees. The mixed coniferous tree (p < 0.001) and coniferous tree (p < 0.05) demonstrate distinctions from Hc mosses in the same context. These findings emphasise the influence of moss species and forest cover on moss collection in the region. Quantitative covariates, EMEP-air, show high statistical significance (p < 0.001) with a positive effect, consistent with the broader model encompassing France. Additionally, covariates Forest_1 and Fruits_10 positively influence the outcome (p < 0.05), while Pasture_5 and Sea_10 exhibit negative effects at the same significance level. Including covariates such as Fruits_10, which are absent in the French model, provides additional insights into the dynamics of the Atlantic zone.

<figure><img src="image_6.png" /><figcaption>**Fig. 2** a Predictions of Cd concentration accumulated by *Hypnum cupressiforme* Hedw. under deciduous cover tree (µg/g) from the selected 9 covariate model. White lines, border of the biogeographical zones. b Standard deviation of the prediction from the selected 9 covariate model</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 3 a Ordinary kriging of Cd concentration in mosses (µg/g). b Standard deviation of ordinary kriging</figcaption></figure>

<table><caption>**Table 6** Leave-one-out cross-validation (LOOCV) with root-mean-square deviation (RMSE) to compare the accuracy of models in France and each biogeographical zone. Atlantic, Continental, Mediterranean, and Alpine</caption><thead><tr><th></th><th>LOOCV RMSE<br>Encompassing France</th><th>LOOCV RMSE<br>Biogeographical zone</th></tr></thead><tbody><tr><td>France</td><td><strong>0.2961672</strong></td><td><strong>0.3620656</strong></td></tr><tr><td>Atlantic zone</td><td><strong>0.2542687</strong></td><td><strong>0.2467906</strong></td></tr><tr><td>Continental zone</td><td><strong>0.3290859</strong></td><td><strong>0.3239615</strong></td></tr><tr><td>Mediterranean zone</td><td><strong>0.3766574</strong></td><td><strong>0.2872906</strong></td></tr><tr><td>Alpine zone</td><td><strong>0.3250785</strong></td><td><strong>0.9370611</strong></td></tr></tbody></table>

In the Continental zone, based on 147 observations, the backward-forward AIC selection retains 12 variables in the model. Only 5 have a significant effect on the Cd concentration. Significant qualitative covariates include moss species and tree cover, which are retained in the model. Quantitative covariates such as Vegetation_5, Pasture_5 and Altitude are significant ($p < 0.05$). Two novel covariates, previously unselected by models encompassing France and the Atlantic region, emerge as notable contributors: Vegetation_5, which exhibits a positive effect and altitude, which manifests a negative impact. Geographical considerations can elucidate the inclusion of these covariates. Areas of higher altitude, such as the Jura and Massif Central, substantiate the introduction of the altitude covariate. The higher we are, the less Cd there is. Conversely, Vegetation_5 is linked to the prevalence of natural pasturelands, heaths and shrublands within the continental zone. Additionally, extensive agricultural landscapes in this region are associated with Pasture_5. These two covariates increase the concentration of Cd.

The selection of these covariates in the Continental zone reflects a logical and context-specific rationale, aligning with this territory's distinctive environmental characteristics and land use patterns. Consequently, the absence of specific covariates, like Sea_10, is expected, as they are unrelated to the specific attributes of this region.

In the Mediterranean zone, based on 49 observations, we applied the backward-forward AIC selection, which retains 20 variables. 11 of them have a significant effect on the Cd concentration. The qualitative covariates corresponding to the tree cover are significant, while moss species do not show significance in this zone. New quantitative covariates are selected and significant in this model: Pasture_1, Vegetation_1, Water_15, Sea_5 and Road_5. These covariates, previously unselected in models across France, the Atlantic and the Continental zones, contribute new information specific to the Mediterranean zone. EMEP-air, Railway_5, Vegetation_5, agricultural_5 and Forest_5 are the same as the previous region and France (Table 4). The positive effect of Railway_5, contrasting with its effect in the Continental zone, is potentially linked to factors such as Mediterranean forest fires (Hubert et al 2002). Moreover, Railway_10 and Road_5, two road-related covariates, demonstrate statistically significant coefficients with a negative impact on Cd concentration, likely due to the extensive railway and road networks in the Mediterranean zone. Several other covariates, including Vegetation_1, Water_15, Sea_5, Vegetation_5, Agricultural_5 and Forest_5, are consistent with the environmental characteristics unique to the Mediterranean zone. The covariate associated with atmospheric conditions, EMEP-air, consistently exhibits a positive effect in this region, as observed in other zones where it is included.
Among the highly significant covariates (*p* < 0.01), Paturelage_1, Vegetation_1, Ferre_5 and Routier_5 have positive coefficients, while Vegetation_5 has a negative coefficient.

In the Alpine zone, our observations are constrained by a limited sample size of 29 against the number of variables (22 variables retained by the AIC selection). However, only three covariates are found to affect Cd concentration significantly. The qualitative covariate of moss species is notably significant. In quantitative covariates, only two show statistical significance: EMEPDep1 and EMEP-air. EMEPDep1 has a negative effect, while EMEP-air has a positive effect across most regions in France, except for the Continental zone, where it has no significant impact. This suggests that EMEP-air contributes to increased Cd concentration across the various regions.

As for encompassing France, we examine the correlation between observed Cd concentrations before and after accounting for covariates in biogeographical zones to assess their impact on spatial correlation in Cd mosses concentrations (response variable). Table 5 provides a summary of these findings.

We are pleased to see that covariates effectively address spatial correlation in Cd concentrations in mosses across three biogeographical zones (Continental, Mediterranean, Alpine).

On the other hand, in the Atlantic zone, the regression model without covariates yields a highly significant Moran I Index (*p* < 0.05) and an I statistic of 0.424. Conversely, the model incorporating 10 variables leads to a less significant *p*-value and an I statistic of 0.011. Notably, including

covariates in this zone substantially account for the spatial correlation of Cd concentrations observed in mosses. This improvement surpasses that observed for the encompassing French model, for which the spatial correlation is only halved. The predictions of our maps will likely be completely independent.

## Prediction map

From each final model selected in each zone, we construct a Cd prediction map for the biogeographical zone (and a map of the standard deviations of the predictions). Figure 4a depicts the composite map of these regional predictions, culminating in a comprehensive depiction of France based on distinct biogeographical models. This is a much more nuanced map of France in each biogeographical zone than the map encompassing France (Fig. 2a) and the kriging (Fig. 3a).

The Cd prediction in the Atlantic zone shows higher Cd concentrations near Lille (1 µg/g Cd) and Paris (0.65 µg/g Cd) regions in northern France. A similar pattern emerges when we compare these results with the model encompassing France (0.65 µg/g) near Lille and Paris. However, the model covering the Atlantic zone offers a more nuanced perspective, incorporating new covariates that enhance the depth of information and our understanding of the observed patterns. By comparing our results visually with the kriging map, we note that areas near Lille and Paris are also present. Our model does indeed find this spatial similarity. The Atlantic zone shows the highest standard deviation

<figure><img src="image_8.png" /><figcaption>**Fig. 4** a Predictions of Cd concentration accumulated by *Hypnum cupressiforme* Hedw. under deciduous cover tree (µg/g) from each model of biogeographical zones. White lines, border of the biogeographical zones.</figcaption></figure>

b Standard deviation of the prediction from each model of the biogeographical zone

<figure><img src="image_9.png" /></figure>
near Paris, corresponding to higher Cd prediction values in mosses. This phenomenon is not found near Lille along the Belgian border, which has more observations (Fig. 1a). The standard deviation map shows greater variability where predictions of high and low concentrations are found.

The Continental zone shows slightly lower Cd concentrations along railway lines, as the covariate is present in the model but not significant (the standard deviation of these predictions is high in these areas, as evidenced in Fig. 4b). We can see from Fig. 3a that the railway zones are absent if we visually compare our results for this zone with the kriging map. However, the covariate coefficient is not significant in the model (Table 3).

In the Mediterranean zone, very high Cd concentrations are predicted near railways within a 5-km buffer due to the significant positive effect of this covariate in the model and the important presence of railways in the region (many predictions in the 2 × 2-km grid have a very high railways covariate value). In addition, we can observe the influence of the negative effect of the Sea_5 covariates along the west coast near the city of Marseille with lower Cd concentration predictions. The highest values modelled by EMEP\_air are also represented near Marseille. EMEP\_air is significant in the model, explaining the high values found in this area. Thus, the Mediterranean zone map offers a new perspective compared with the other maps (Fig. 2a and Fig. 3a), highlighting the possibility of high Cd concentrations in mosses along railway lines. However, the standard deviation of these predictions is high.

The Alpine zone highlights a significant impact of EMEP\_air, which is present in the regression model with a positive effect. Let us compare this result visually with the map of encompassing France. We note that the impact of the EMEP\_air covariate is absent from the predictions constructed on the encompassing model for France in this zone, whereas the kriging map reveals some of this effect.

Table 6 presents the LOOCV scores of the biogeographical model calculated on a zone-by-zone basis and France using the different models. Employing separate models for each zone for the Atlantic, Continental and Mediterranean zones yields superior LOOCV scores compared to using a single model for all of France. However, in the Alpine zone, the LOOCV score for the model encompassing all of France (0.328) is lower than the zone-specific LOOCV score (0.937). This discrepancy is observed since the model has nearly as many covariates as observations, and removing an observation to calculate an LOOCV criterion probably makes it unreliable. Consequently, the overall LOOCV score for France, calculated using different models for each zone (0.362), surpasses that obtained with a single model for the entire country (0.2961). Nevertheless, the LOOCV score for France using the zone-by-zone approach remains lower than that obtained through the kriging method (0.491).

# Discussion

## Covariates selected and comparison with literature

Quantitative covariates such as urban, agricultural, sea or forest land use were categorised from the CLC datasets (Nickel et al. 2017). We used numerous transformations based on the CLC (Table 1) to leverage the most relevant covariates strengthening their link to Cd concentrations in mosses. For instance, while the CLC subclass ‘grove’ was denoted by the 222 and 223 codes, we introduced another category labelled ‘fruits’ including these two codes as well as the 221 code which corresponded to the ‘wine-growing’ label.

The present study sought to analyse relevant covariates such as moss species, urban land uses, altitude and Cd air concentration EMEP modelled data, which were assessed in prior research by Schröder et al. (2010), Nickel et al. (2017), Harmens et al. (2014) and Holy et al. (2009). Notably, although EMEP encompasses Cd deposition and air concentrations, previous studies predominantly focused on EMEP Cd deposition. In this analysis, we included deposition and air concentration metrics, revealing that air concentration emerged as the preferred covariate based on its significant *p*-values, except for the Continental zone. A study considering all of Europe (Holy et al. 2009) has modelled these two EMEP covariates. In this study, the covariate selection retained the EMEP deposition covariate as significant. Another relevant study focusing solely on Germany (Schröder et al. 2019) also identified EMEP deposition as a significant covariate. This result could be due to multiple geographical and biological differences in the predefined territories of Europe, Germany and France. Our study globally examined a wider range of variables (22 in the Alpine zone to 29 in the Atlantic zone) compared to previous studies by Schröder et al. (2010) and Nickel et al. (2017), investigating 11 and 9 variables, respectively. In line with prior research, the EMEP coefficient emerged as a significant predictor in our analysis. However, we noted spatial dependence in the EMEP coefficient, exhibiting varying values across different regions of France. The observed spatial variability underscores the importance of conducting further investigations into the complex relationships between EMEP and Cd concentrations in moss. Moreover, Cd air concentrations were distributed across three French biogeographical zones, while Cd deposition was evidenced only in the Alpine and Continental zones.

Between 2010 and 2021, major sources of Cd were identified as originating from industries, transportation, human activities, agricultural practices and waste, as reported in Citepa (2023). Our model established a link
between Cd concentrations in moss in France and the fol-
lowing selected covariates: urban, industrial and agricul-
tural land uses. Parallel studies conducted by Lequy et al.
(2022) and Schröder et al. (2017), Schröder et al. (2010))
also highlighted the relevance of these covariates, particu-
larly the urban and agricultural influences. Furthermore,
the divergence from the biogeographical model was appar-
ent in our study across distinct Atlantic, Continental, Med-
iterranean and Alpine zones. Consequently, specific pat-
terns emerged when examining the biogeographical zones.
The Atlantic zone showed associations with relevant fac-
tors such as fruit (relevant to wine-growing), which had
a positive effect, and pasturelands, which has a negative
effect. The Mediterranean zone exhibited stronger correla-
tions with transportation networks (railway and road infra-
structure). Moreover, qualitative covariates such as moss
species, tree cover and biogeographical zone consistently
influence the French model. These outcomes resonate with
a prior study (Lequy et al. 2017a, b), highlighting the con-
sistent selection of similar moss species and quantitative
covariates like Cd air concentration, industrial land uses
and urban land uses.

Selection criteria

Diverse modelling approaches, such as linear regression,
linear mixed models and decision trees, were employed to
assess these covariates. Several alternative methodologies
were considered, such as classification and regression trees
(CART) (Holy et al. 2009) and linear regression with R²
selection (Harmens et al. 2011). The difference lies in the
selection constraints. For instance, while R² selection signifi-
cantly improved model goodness of fit, it did not penalise the
model size (i.e. more covariates led to better goodness of fit).
Our selection process followed a three-step methodology:
starting with Spearman correlation (Schröder et al. 2019),
followed by expert input to carefully limit the number of
variables while retaining relevant information for our model.
Subsequently, we employed the AIC (Lequy et al. 2017a, b)
using a backward-forward selection technique, ensuring an
unbiased model selection process. We reduced dimensional-
ity and minimised the AIC while retaining non-significant
covariates to enhance model goodness of fit. Given the
large number of covariates at the start of the AIC selection
(around thirty), an alternative might have been to leverage
a ridge or lasso method. The ridge or lasso method is used
in the context of a high number of variables relative to the
size of the number of observations (Tibshirani 1996). We
therefore applied these methods in our case study. However,
these approaches were either too sparse (e.g. lasso did not
retain any covariates in the biogeographical zone) or too
inclusive (e.g. ridge does not eliminate any covariates in
the biogeographical zone). Thus, we found the AIC to be a

suitable compromise between model parsimony and good-
ness of fit, consistent with existing literature.

Some studies on mosses, such as those conducted by Holy
et al. (2009) and Schröder et al. (2019), often overlooked the
spatial dimension in their models when selecting and apply-
ing covariates. Our investigation identified spatial correla-
tion in Cd concentrations across France, a factor discarded
during covariate selection in previous studies. To address
this issue, we initially partitioned France into biogeographi-
cal zones, which mitigated the problem. As a result, residual
correlation in 3 out of 4 models was negligible after con-
sidering the covariates. A spatial statistical model could be
applied to fully assess the spatial pattern of the data (Dor-
mann et al. 2007).

Relevance of biogeographical zone division

As mentioned previously, the partition of France in the bio-
geographical zone has already demonstrated the benefits of
introducing relevant covariates in each zone to elucidate
the spatial correlation over France. Another advantage was
shown when evaluating the leave-one-out cross-validation
(LOOCV) criterion to assess the predictive performance of
each zone. In three out of four zones, these criteria dem-
onstrated improved performance compared to the analysis
encompassing the entirety of France. Although this zon-
ing strategy enhanced model performance, its reliance on
a limited number of observations raised concerns regarding
the reliability of measures such as LOOCV. This limitation
was particularly evidenced in the Alpine zone, which com-
prised 29 observations and 22 variables selected via AIC
backward-forward selection. An alternative partitioning
of West Europe, known as the Ecological Land Classifica-
tion of Europe (ELCE) (Hornsmann et al. 2008), provided
a division into 40 distinct regions, which gave us a finer
breakdown of France (more than 4 zones). Such a delimita-
tion may prove to be beneficial in scenarios where a larger
volume of observations is available for analysis.

Benefits and drawbacks of multiple linear
regression models to make prediction maps

Applying linear models for predicting Cd concentration
in mosses produced cartographic representations with
enhanced contrast and more comprehensive information
content than kriging maps. Each significant coefficient in
our model played a crucial role in shaping the structure of
the resulting map. Interpretation of these maps was facili-
tated by examining the estimated coefficients associated
with the relevant covariates (Schröder et al. 2019). For
instance, the Mediterranean zone became discernible in
our model which may not be the case in kriging maps. This
is due to the inclusion of significant covariates such as
EMEP\_air or Sea\_5 in our model. The presence of specific covariates, such as Sea\_5 reduced Cd concentration in the Mediterranean zone. This phenomenon can be attributed to the sea salt effects (Renaudin et al. 2018) but also to the direction of the prevailing winds in this region. In the French metropolis, the winds come mainly from the land, even though the sites are close to the sea, while in Corsica, they are arriving from the sea. Furthermore, establishing a connection between the modelisation of atmospheric Cd concentration (EMEP\_air) and the concentrations observed in mosses added valuable insight. It has to be noted that the validity of these predictive maps relied on the validation of moss measurements within these prediction zones. However, regression models exhibited limitations (Zuur et al. 2009), especially when dealing with a limited number of observations, which can decrease their reliability. These models may face challenges in accommodating outliers, leading to the generation of extreme values. Additionally, predicting a zone posed specific challenges due to variations in environmental conditions. The transition from collecting mosses in areas under tree cover, away from pollution sources, to predicting in urban zones introduces a stark environmental contrast. This disparity can challenge achieving coherent predictions, potentially leading to extreme values. To address this issue, we considered transforming the response variable, imposing constraints that aligned with the characteristics of the collected mosses. This transformation helped to ensure that the predictions remained within a meaningful and realistic range of values, accounting for the distinct conditions encountered in different zones. Employing a complementary log-log transformation was a viable strategy to circumvent these issues. Despite confining the predictive values within the Cd concentration range, this transformation limited our ability to derive mean-based maps directly. As a result, we generated maps based on medians. To assess the variability of our predictive models, we also examined the standard deviation derived from the model predictions. Analysing the standard deviation offered valuable insights to quantify the model's ability to capture and explain the observed data.

Furthermore, we generated a spatial grid with a 2×2-km resolution to calculate our predictions based on a multiple regression complementary log-log model. The model was constructed using concentrations collected from mosses situated under forest cover, away from local sources of pollution across France. Our predictions spanned 139,673 pixels on the grid throughout the country, with a varying number of pixels in different biogeographical zones. Notably, this grid enabled predictions in urban zones close to pollution sources. The choice of this fine-scale grid allowed us to produce detailed maps that convey

richer information due to the inclusion of covariates, as opposed to traditional kriging methods.

The LOOCV criterion was a valuable tool for looking at the prediction performance. It involved systematically removing one observation from the model and then predicting it using the remaining data. This process assessed the ability of a model to predict new data points. However, it did not necessarily reflect its performance within the observed data space. Moreover, caution is needed when extrapolating beyond the observed data, especially in cases where the prediction grid includes geographic areas significantly different from the observed data, such as urban regions. Further measurements within the territory can provide full validation of our predictions.

## Conclusion

Using 445 sites in a multiple linear regression model to predict and visualise Cd background levels in mainland France provided valuable insights into the factors influencing Cd concentrations across the French territory. The resulting models also allowed for the creation of informative prediction maps on a fine scale of the territory. It highlighted the relationships between Cd concentrations in mosses and various anthropogenic environmental factors, such as railroads, and the concentrations of airborne Cd accumulated by mosses. Our method produced more detailed maps compared to those generated by ordinary kriging. Adding a spatial effect could increase its accuracy in order to improve this prediction and mapping model.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s11356-025-35985-5.

**Acknowledgements** The authors wish to acknowledge all the collectors who contributed to the BRAMM survey and those who took care of the samples in the laboratory. The authors would like to thank the staff of the USRAVE laboratory (INRAE- Centre de Bordeaux), accredited by The French Committee for Accreditation, for their analyses of Cd accumulated by mosses. The authors thank all field and labworkers involved in the RMQS and data managers. The “Groupement d’intérêt scientifique sol” funded soil sampling and classical soil analysis within the French National soil monitoring network framework. The authors are also thankful to Ilia Ilyin of MSC-East of EMEP for providing the modelled atmospheric deposition and concentration of Cd.

**Author contribution** All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Jérémy Lamouroux, Caroline Meyer, Sébastien Leblond and Isabelle Albert. The first draft of the manuscript was written by Jérémy Lamouroux, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

**Funding** Open access funding provided by Université Paris-Saclay. The French Agency for Ecological Transition (conventions ADEME-MNHN N°1562c0027 and 1666C0022) supported Mosses collection and metal concentration analysis. Jérémy Lamouroux received a PhD
grant from the National Research Institute for Agriculture, Food, and the Environment (INRAE) in France and The French National Museum of Natural History (MNHN).

**Data availability** Data on Cd concentrations in mosses are publicly available on request from Sébastien Leblond (sebastien.leblond@mnhn.fr) or Caroline Meyer (caroline.meyer@mnhn.fr). All models used for this work are publicly available. The localisation of the data can be accessed via the following dataset on the Global Biodiversity Information Facility (GBIF): https://www.gbif.org/dataset/166cbacd-528a-4a84-b1db-0177b50ac946. Information regarding the Cd measurements is available at: https://bramm.mnhn.fr/cadmium-cd/.

## Declarations

**Ethical approval** Not applicable for that specific section.

**Consent to participate** Not applicable for that specific section.

**Consent for publication** Not applicable for that specific section.

**Competing interests** The authors declare no competing interests.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Andujar P, Bensefa-Colas L, Descatha A (2010) Acute and chronic cadmium poisoning. Rev Med Interne 31:107-115. https://doi.org/10.1016/j.revmed.2009.02.029

Banerjee S, Carlin BP, Gelfand AE (2003) Hierarchical modeling and analysis for spatial data. Chapman and Hall/CRC. https://doi.org/10.1201/9780203487808

Bates JW (1992) Mineral nutrient acquisition and retention by bryophytes. Jbryology 17:223-240. https://doi.org/10.1179/jbr.1992.17.2.223

Berg T, Steinnes E et al (1997) Use of mosses (Hylocomium splendens and Pleurozium schreberi) as biomonitors of heavy metal deposition: from relative to absolute deposition values. Environ Pollut. https://doi.org/10.1016/S0269-7491(97)00103-6

Boquete MT, Ares A, Fernández JA, Aboal JR (2020) Matching times: trying to improve the correlation between heavy metal levels in mosses and bulk deposition. Sci Total Environ. https://doi.org/10.1016/j.scitotenv.2020.136955

Butt A, Qurat-ul-Ain RK et al (2018) Bioaccumulation of cadmium, lead, and zinc in agriculture-based insect food chains. Environ Monit Assess. https://doi.org/10.1007/s10661-018-7051-2

Citepa (2023) Gaz à effet de serre et polluants atmosphériques. Bilan des émissions en France de 1990 à 2022. Rapport Secten éd. 2023 (in French). Available on 21 March 2024 at https://www.citepa.org/wp-content/uploads/publications/secten/2023/Citepa_Secten_ed2023_v1.pdf

Cressie N (1990) The origins of kriging. Math Geol 22:239-252. https://doi.org/10.1007/BF00889887

Dormann CF et al (2007) Methods to account for spatial autocorrelation in the analysis of species distributional data: a review. https://doi.org/10.1111/j.2007.0906-7590.05171.x

Gurjar BR, Molina LT, Ojha CSP (2010) Air pollution: health and environmental impacts (1st ed.). CRC Press. https://doi.org/10.1201/EBK1439809624

Gusev A, Iliyn I et al (2009) Modelling of heavy metals and persistent organic pollutants: new developments (Technical Report No. 1/2009). EMEP/MSC-E, Moscow, Russia. Available on 21 March 2024 at https://www.msceast.org

Harmens H et al (2011) Nitrogen concentrations in mosses indicate the spatial distribution of atmospheric nitrogen deposition in Europe. Environ Pollut. https://doi.org/10.1016/j.envpol.2011.04.041

Harmens H et al (2012) Country-specific correlations across Europe between modelled atmospheric cadmium and lead deposition and concentrations in mosses. Environ Pollut. https://doi.org/10.1016/j.envpol.2012.02.013

Harmens H et al (2014) Relationship between site-specific nitrogen concentrations in mosses and measured wet bulk atmospheric nitrogen deposition across Europe. Environ Pollut. https://doi.org/10.1016/j.envpol.2014.07.016

Holy M, Leblond S, Pesch R et al (2009) Assessing spatial patterns of metal bioaccumulation in French mosses by means of an exposure index. Environ Sci Pollut Res. https://doi.org/10.1007/s11356-009-0146-0

Hornsmann I et al (2008) Calculation of an Ecological Land Classification of Europe (ELCE) and its application for optimising environmental monitoring networks. In: Geospatial Cross-roads@ GL_Forum'08: Proceedings of the Geoinformatics Forum Salzburg. Wichmann, Heidelberg: 140-151.

Hubert T et al (2002) Plans de prévention des risques naturels (PPR) Risques d'incendies de forêt. Ministère de l'écologie, France. Available on 21 March 2024 at https://www.ecologie.gouv.fr/sites/default/files/publications/PPR feux de foret complet 0.pdf

Ilyin I et al (2013) Heavy metal transboundary pollution of the environment (No. Status Report 2/2013). EMEP/MSC-E/CCC/CEIP. Available on 21 March 2024 at https://msceast.org/reports/2_2010.pdf

Lequy E et al (2017a) 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. Environ Pollut. https://doi.org/10.1016/j.envpol.2016.10.065

Lequy E et al (2017b) Spatial analysis of trace elements in a moss bio-monitoring data over France by accounting for source, protocol and environmental parameters. Sci Total Environ. https://doi.org/10.1016/j.scitotenv.2017.02.240

Lequy E et al (2019) Long-term exposure to atmospheric metals assessed by mosses and mortality in France. Environ Int. https://doi.org/10.1016/j.envint.2019.05.004

Lequy E et al (2022) Influence of exposure assessment methods on associations between long-term exposure to outdoor fine particulate matter and risk of cancer in the French cohort Gazel. Sci Total Environ. https://doi.org/10.1016/j.scitotenv.2022.153098

Lequy E et al (2023) Long-term exposure to airborne metals and risk of cancer in the French cohort Gazel. Environ Int. https://doi.org/10.1016/j.envint.2023.107999

Nickel S et al (2017) Modelling and mapping heavy metal and nitrogen concentrations in moss in 2010 throughout Europe by applying random forests models. Atmos Environ. https://doi.org/10.1016/j.atmosenv.2017.02.032

Nordberg GF et al (2007) Handbook on the toxicology of metals. Elsevier. https://doi.org/10.1016/B978-0-12-369413-3.X5052-6
Renaudin M et al (2018) The coastal environment affects lead and sodium uptake by the moss Hypnum cupressiforme used as an air pollution biomonitor. Chemosphere. https://doi.org/10.1016/j.chemosphere.2017.11.045

Roekaerts M (2012) Future of the European diploma for protected areas state of network, analysis of the different types of habitats and biogeographic regions already represented and proposals for developing the network into the future, convention on the conservation of European wildlife and natural habitats. Available on 21 March 2024 at https://rm.coe.int/1680746283

Schober P et al (2018) Correlation coefficients: appropriate use and interpretation. Anesth Analg. https://doi.org/10.1213/ANE.00000 00000002864

Schröder W et al (2010) Are cadmium, lead and mercury concentrations in mosses across Europe primarily determined by atmospheric deposition of these metals? J Soils Sediments. https://doi.org/10.1007/s11368-010-0254-y

Schröder W et al (2014) Mapping correlations between nitrogen concentrations in atmospheric deposition and mosses for natural landscapes in Europe. Ecol Ind. https://doi.org/10.1016/j.ecolind.2013.09.013

Schröder W et al (2017) Bioindication and modelling of atmospheric deposition in forests enable exposure and effect monitoring at high spatial density across scales. Ann for Sci. https://doi.org/10.1007/s13595-017-0621-6

Schröder W et al (2019) Moss species-specific accumulation of atmospheric deposition? Environ Sci Europe. https://doi.org/10.1186/s12302-019-0262-5

Shim H, Bonifay W, Wiedermann W (2023) Parsimonious asymmetric item response theory modeling with the complementary log-log link. Behav Res Methods. https://doi.org/10.3758/s13428-022-01824-5

Stone M (1974) Cross-validatory choice and assessment of statistical predictions. J Roy Stat Soc: Ser B (Methodol). https://doi.org/10.1111/j.2517-6161.1974.tb00994.x

Taeprayoon P et al (2023) Potentially toxic element accumulation of bryophyte taxa in contaminated soils at Tak Province, Thailand. Ecol Indic. https://doi.org/10.1016/j.ecolind.2023.109971

Tibshirani R (1996) Regression shrinkage and selection via the lasso. J Roy Stat Soc: Ser B (Methodol). https://doi.org/10.1111/j.2517-6161.1996.tb02080.x

Travnikov O, Ilyin I (2005) Regional model MSCE-HM of heavy metal transboundary air pollution in Europe. EMEP/MSC-F Technical Report 6/2005, p 59. Available on 21 March 2024 at https://d1wqtxts1xzle7.cloudfront.net/83539425/Regional_Model_MSCE-HM_of_Heavy_Metal-Tr20220409-3729-q71uko.pdf?16494

89187=&response-contentdisposition=inline%3B+filename%3DRegional_Model_MSCE_HM_of_Heavy_Metal-Tr.pdf&Expires=1738232268&Signature=Vuk0cXKZLfKb87f4t7-h0AtKSAf4vLlTH1bS6~9eE3EIfiz6M39ANFmPgGG~2f1ml8bu0Jqi1lj bUz3P5bToYTraQaybe3IPvsurf03uCp5jKHh9q476eMnoffgcV47dzfzPS0sslhxbvzqjVf1TnHFjPhKdqWKQYRUNNOTDN1CqAw-C3q4povJvqXjSuLmp~lG6wjRHzphJyQ1YIUYFG~Gngio3vCix6MJSuUgAq2-08YV~MROqL0sDLshRDkaOloP6Y7SPzYyL6~UdCytDZ4dTБ5vzY2KVMnimsIg6Cznv2VZnh2P7Q4dwCvRNPOZs04hp3-WL9EP3UXwdUh7tg__&Key-Pair-Id=APKAJLOHF5GGSLRBV4ZA

Tyagi K et al (2022) Artificial intelligence and machine learning for EDGE computing. Academic Press, Chapter 4 - Regression analysis, pp 53–63, https://doi.org/10.1016/B978-0-12-824054-0.00007-1

Vieille B et al (2021) Are Grimmia mosses good biomonitors for urban atmospheric metallic pollution? Preliminary evidence from a French case study on cadmium. Atmosphere. https://doi.org/10.3390/atmos1204091

Wilcox R et al (2012) Introduction to robust estimation and hypothesis testing (third edition). https://doi.org/10.1016/C2010-0-67044-1

Zhang H, Reynolds M (2019) Cadmium exposure in living organisms: a short review. Sci Total Environ. https://doi.org/10.1016/j.scitotenv.2019.04.395

Zuur AF, Ieno EN, Walker NJ, Saveliev AA, Smith GM (2009) Limitations of linear regression applied on ecological data. In: Mixed effects models and extensions in ecology with R. Statistics for Biology and Health. https://doi.org/10.1007/978-0-387-87458-6_2

Zwart D et al (2010) Loss of species due to cadmium and lead depositions in Europe. Progress in the Modelling of Critical Thresholds and Dynamic Modelling, including Impacts on Vegetation in Europe, 101. Available on 21 March 2024 at https://rivm.openg repository.com/bitstream/handle/10029/258798/680359001.pdf?sequence=3#page=103

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

The submitted work is original and has not been published elsewhere in any form or language (partially or in full). The results presented are without fabrication, falsification or inappropriate data manipulation (including image-based manipulation).