Trace metals from historical mining sites and past metallurgical activity remain bioavailable to wildlife today

Throughout history, ancient human societies exploited mineral resources all over the world, even in areas that are now protected and considered to be relatively pristine. Here, we show that past mining still has an impact on wildlife in some French protected areas. We measured cadmium, copper, lead, and zinc concentrations in topsoils and wood mouse kidneys from sites located in the Cévennes and the Morvan. The maximum levels of metals in these topsoils are one or two orders of magnitude greater than their commonly reported mean values in European topsoils. The transfer to biota was effective, as the lead concentration (and to a lesser extent, cadmium) in wood mouse kidneys increased with soil concentration, unlike copper and zinc, providing direct evidence that lead emitted in the environment several centuries ago is still bioavailable to free-ranging mammals. The negative correlation between kidney lead concentration and animal body condition suggests that historical mining activity may continue to play a role in the complex relationships between trace metal pollution and body indices. Ancient mining sites could therefore be used to assess the long-term fate of trace metals in soils and the subsequent risks to human health and the environment.


Results
Trace metal concentrations in topsoils. The TM concentrations in soils ranged from less than 0.5 mg·kg −1 to 54.2 mg·kg −1 for Cd, 11 mg·kg −1 to 212 mg·kg −1 for Cu, 85 mg·kg −1 to 8410 mg·kg −1 for Pb, and  (Table 1 and Supplementary Table S1). Maximum levels for all the four TMs studied were found in the contaminated sites of the Morvan region. The spatial distribution of TMs in soils shows that higher concentrations were found in the mining and metallurgical sites, whatever the region considered ( Fig. 2a, Fig. 3a, Supplementary Fig. S1 for the Morvan and Supplementary Fig. S2 for the Cévennes). In the Morvan, the three sites differed significantly (Kruskal-Wallis test, all p < 0.05) in terms of Cd, Cu, Pb and Zn contents in topsoils ( Fig. 2a and Fig. 3a), with all elements following the pattern M0 < M1 < M2 (Steel-Dwass pairwise comparisons). In the Cévennes, the three sites differed for Cu and Zn contents with the pattern C2 < C0 < C1 (Fig. 2a). The Cd concentrations were similar in C0 and C2, both were lower than in C1. The Pb concentrations were similar for the three study sites (Kruskal-Wallis test, p > 0.05) but C2 exhibited greater spatial heterogeneity, ranging from 31 mg·kg −1 to 4810 mg·kg −1 ( Fig. 3a and Supplementary Table S1).
Trace metal concentrations in kidneys. The Cu concentrations in wood mouse kidneys ranged from 10.37 μg·g −1 to 22 μg·g −1 , and the Zn concentrations ranged from 44.5 μg·g −1 to 160 μg·g −1 (Supplementary  Table S2). Except for Cu in the Cévennes with C1 < C0 − C2 (Tukey's HSD test after an analysis of variance (ANOVA) with p < 0.05), essential element concentrations did not differ between the study sites (p > 0.05), suggesting physiological (homeostatic) regulation (Fig. 2b). Concentrations of the non-essential elements in the kidneys ranged from 0.05 μg·g −1 to 38 μg·g −1 for Cd and from 0.05 μg·g −1 to 19 μg·g −1 for Pb. Maximum concentrations were found in the contaminated sites of the Cévennes region (C1 for the Cd and C2 for Pb). The Cd concentrations followed these patterns: M0 < M1 − M2 and C0 − C2 < C1. The Pb concentrations followed the same pattern in both parks, showing no differences between sites 0 and 1, but with both values lower than for site 2 (Fig. 3b).

Trace metal concentrations in wood mice in relation to biological and environmental parameters.
Multivariate linear models were used to investigate the relationship between TM concentrations in kidneys and explanatory variables (site, TMs in soils, sex, and mass), for each TM separately. Both Cd and Zn concentrations in wood mouse kidneys were best explained by models combining study sites and body mass ( Table 2, for  description of best-fit models and Supplementary Table S3 for model parameters). The same pattern was observed for Cu concentrations in wood mouse kidneys, with sex as an additional factor (Table 2 and Supplementary  Table S3). The Pb concentrations varied between sites as indicated above and increased with Pb concentrations in soil (Table 2, Fig. 4a and Supplementary Table S3). The Cd concentrations in wood mouse kidneys increased slightly with mass, while Cu and Zn concentrations decreased slightly (Table 2 and Supplementary Table S3).

Toxic effects assessed by body condition, somatic indices, and fluctuating asymmetry.
Relationships between body condition and somatic indices were investigated with multivariate linear models using study sites, the four TMs in wood mouse kidneys, and sex as explanatory variables. Body condition as assessed by scaled mass index (SMI) was best explained by a model combining study sites, Cd and Pb concentrations in kidneys and the interaction between Cu or Zn concentrations in kidneys and sex ( Table 3). The SMIs varied according to the sites, increased with Cd concentrations and were negatively related to Pb concentrations in kidneys ( Fig. 4b and Supplementary Table S4). The SMIs were negatively influenced by the interaction between Cu concentrations and sex, and positively influenced by the interaction between Zn concentrations and sex. Somatic index data (scaled liver index -SLI, and scaled kidneys index -SKI) were best fit by models including study sites, sex and the interaction between Cu or Zn concentrations in kidneys and sex (Table 3). Like SMI, both SKI and SLI were negatively related to the Cu × sex interaction and positively related to the Zn × sex interaction (Supplementary Table S4). Concerning FA, preliminary tests were performed for all traits measured (here, length Cd (mg · kg −1 ) Cu (mg · kg −1 ) Pb (mg · kg −1 ) Zn (mg · kg −1 ) and width of lower molars) as recommended by Palmer 29 . These tests did not suggest that directional asymmetry (DA), antisymmetry (AS) or relationships between asymmetry and trait size could significantly bias FA estimates (Supplementary Table S5). No significant correlations were found between absolute asymmetry distribution and TM concentrations in kidneys (Supplementary Table S6). When compared to measurement error, FA10, a parameter used for FA assessment, was always significant (Supplementary Table S7), but no clear relationship was found between levels of developmental instability assessed by FA and sites for either park ( Supplementary Fig. S4).

Discussion
Even centuries after mining and metallurgical activities have ceased, TM concentrations in soils surrounding such sites still reach high levels. There is no European consensus on threshold concentrations of metals in soils, but the Dutch government proposed a soil classification scheme in which intervention values are defined. These strict values identify serious contamination of soils (12,190,530, and 720 mg·kg −1 dry matter for Cd, Cu, Pb, and Zn, respectively) and indicate that remediation is necessary 31 . In the present study, 75 out of 261 soil samples (29%) exceed the Pb intervention value, reaching 59% in the highly contaminated Morvan site, M2 (Fig. 3a). For both Cd and Cu, only one sample out of 261 exceeds the Dutch intervention values, while Zn exhibits an intermediate pattern, with 10% of the soils exceeding the threshold (Fig. 2a). Among the four metals studied and according to this classification, Pb concentrations represent the most important risk, which is probably linked to the nature of the ore exploited (galena) and to the low mobility of Pb in the environment 32 . For Cd, while only one soil exceeds the Dutch intervention threshold, this intervention limit is high (12 mg·kg −1 ) compared to the average concentration in European surface soils (0.28 mg·kg −1 , Table 1) 32 . In France, Cd concentrations in sewage sludge-amended agricultural soils must be lower than 2 mg·kg −1 ( Table 1). As Cd is toxic to organisms at low doses, it should also   Fig. S1 and Fig. S2), complicating risk assessment. Such heterogeneity is also observed in modern mining or smelting sites 33,34 and can be explained by the spatial distribution of exploitation facilities in relation to the surrounding habitat (e.g., interception of metals emitted in the air by the canopy 35 ).
Metal toxicity depends not only on the concentration of a substance in a medium but also on its bioavailability, a complex combination of the physico-chemical characteristics of the pollutant in abiotic compartments and the biological characteristics of organisms 36 . The first component, known as "environmental availability", represents the physico-chemically driven desorption processes that determine the mobile proportion of the total metal   Table 3. Summary of best-fit models for body condition and somatic indices. Models relating body condition and somatic indices to biological and environmental parameters and trace metal concentrations in wood mouse kidneys. 0 < *** < 0.001 < ** < 0.01 < * < 0.05 < . < 0.1. concentration in a soil. Among the numerous parameters that determine this availability, time generally lead to immobilisation of metals, a process named "ageing" 22,23 . The second component, known as "environmental bioavailability", represents physiologically driven uptake processes that occur when an organism and a pollutant co-occur in time and space (namely, the "exposure" of an organism). In the present work, our aim was to examine whether these high, spatially heterogeneous concentrations of metals were bioavailable to organisms by measuring TM levels in wood mouse kidneys. For the two essential elements, concentrations did not differ among sites, suggesting their efficient (homeostatic) regulation 37 . Concentrations measured here are similar to values found in other studies on wood mice 38,39 . However, the non-essential elements Cd and Pb showed marked differences between sites, with higher concentrations in the mining and smelting areas, showing that these metals are still bioavailable to wildlife. Concentrations of Pb in wood mouse kidneys measured in the present study are of the same order of magnitude as those measured in animals sampled around a non-ferrous smelter still in activity in Antwerp (Belgium) 38 . They are, however, below the concentrations observed in wood mice in the surroundings of Metaleurop Nord, a Pb and Zn smelter in activity from 1894 to 2003 in northern France 28,40 . Direct comparisons between all these values must be undertaken with caution because the sites differ in terms of metal concentrations in soils, soil properties, and/or exposure (diet). Biota-to-soil accumulation factors (BSAFs) 41 40 . This lower value suggests less Pb transfer in the present study, which may indicate lower availability of this metal in soils affected by ancient mining contamination, as suggested by Camizuli et al. 42 .
Once a pollutant has been taken up by an organism, toxic effects will occur only if various defence mechanisms are overcome. In this study, Cd and Pb concentrations in wood mouse kidneys are below the Lowest Observed Adverse Effect Levels (LOAELs), as defined by Shore and Douben 43,44 (Fig. 3b), suggesting that toxic effects are unlikely to occur. However, we found a significant negative relationship between SMI and kidney Pb concentrations. The SMI is a measure of body condition that is often defined as a measure of the energetic (or nutritional) state of an animal 25 . Even if the calculation and interpretation of such indices are still much debated 24,27 , these indices are assumed to be related to fitness. Here, we used the SMI, which has recently been shown to be a better indicator of the relative size of energy reserves than condition indices based on ordinary least squares residuals (see Peig and Green 25 , for details). Although several studies have shown a decline in the body condition of small mammals from polluted sites compared to controls 45,46 , confounding factors like food availability, habitat quality or other chemical elements may contribute to the complex relationships that are observed between pollution and body indices 28 . Therefore the negative relationship between SMI and Pb concentrations, observed in this study, cannot simply be interpreted as implying a direct causal relationship. Other relationships that remain complex to interpret are the positive correlation between SMI and Cd concentrations in kidneys, and the interactions between essential element concentrations (Cu and Zn) and sex. These complex relationships between the body condition of free-ranging vertebrates and both essential and non-essential elements clearly require further investigation. Somatic indices (the relative size of internal organs), which may reveal oedema in individuals exposed to toxic compounds 46 , did not exhibit any clear relationship with metals. Concerning FA, relationships between the degree of developmental instability of populations and the level of the environmental and/or genetic stress to which they were subjected have already been demonstrated 10,47 . In this study, FA10 was detected, but no relationship was found with the sampling site. This result is not in agreement with a study on aquatic ecosystems in the Cévennes National Park, which showed that wild trout were affected by increasing developmental instability in relation to mining contamination 10 .
Taken together these results show that several centuries after mining and smelting activities have ceased, metals are still bioavailable to wildlife, with Pb, and to a lesser extent Cd, concentrations increasing in wood mouse kidneys in relation to soil concentrations. Further studies should be undertaken to determine the precise TM transfer mechanisms that occur in our study sites from the environment to animals (ingestion of soil, animal and plant materials, inhalation of contaminated dust/soil particles, and/or direct transfer through dermal contact). The BSAFs suggest, however, that bioavailability might be lower in soils affected by ancient mining than in soils that have been more recently contaminated. Higher concentrations of Pb in the kidneys of individuals from the most contaminated sites and the negative relationship between these concentrations and SMI raise the issue of the present-day consequences of past anthropogenic activities on wildlife. Specific biomarkers of exposure such as the induction of metallothioneins, a protein involved in the homeostasis of essential elements and in the regulation of non-essential ones 40 , could be envisaged in future studies. Biomarkers of toxic effects, for instance related to the oxidative stress that exposure to trace metals may generate 48 or to histological pathologies 49 , would also provide further insight into the possible toxic effects that may occur in wildlife. Nature Parks are now protected areas and are considered to be relatively pristine, but in the past, they often were the setting for economic and industrial activities. Ancient industrial activities might sometimes have vanished from collective memory but may still represent a risk that deserves investigation. plasma atomic emission spectroscopy (ICP-AES) after pseudo-total aqua regia digestion at Actlabs (Ontario, Canada). This chemical method was chosen because we targeted anthropogenic pollution that is easily extractable, and thus bioavailable to wildlife. Analytical quality control was doubly checked (i) by Actlabs measuring 18 replicates, 8 blanks and several certified reference materials (CRMs) and (ii) by inserting another 25 duplicates and JSD-1, JSD-2, BCSS-1 and PACS-1 CRMs (stream, estuarine, and harbour sediments, respectively) as blind samples. The Actlabs protocol set the limits of detection (LODs) at 0.5 mg·kg −1 for Cd, 1 mg·kg −1 for Cu, 2 mg·kg −1 for Pb and 2 mg·kg −1 for Zn (see Supplementary Tables S8 and S9 for details).

Methods
Small mammal sampling and analysis. Small mammal sampling was conducted on 10 plots of 100 m × 100 m at each site ( Supplementary Figs S1 and S2). These plots were randomly chosen for the two reference sites (M0 and C0) and at the location of anthropogenic activities for the contaminated sites. Wood mice were trapped between mid-September and mid-October 2010. Sampling authorisations were obtained from the DREAL Bourgogne (French regional territory agency in charge) and from the Cévennes National Park. For each of the 10 selected plots, a line of 25 baited traps was set with alternating INRA (door-) and snap-traps, spaced 4 m apart. An extra chamber was added to the classical INRA trap to increase the survival of trapped animals. Traps were set for 3 to 5 consecutive days to ensure an adequate number of samples and were checked and rebaited each morning. The wood mice caught alive were immediately sacrificed by cervical dislocation in accordance with relevant guidelines and regulations 50,51 and frozen as soon as possible after their capture. They were stored at −20 °C until dissection in the laboratory. All small mammals captured were determined at the species level by molecular analysis (sequencing of the cytochrome b gene) of a tissue sample performed by a service provider (ADNid laboratory). Body wet mass was used as an estimator of age 45,52,53 . Combined with reproductive status, three age categories were constructed: juvenile (J), subadult (SA), and adult (A), as in Peig and Green 26 (Supplementary Table S10). For all specimens body mass was measured to the nearest 0.01 g and body length to the nearest 0.01 mm. Livers and kidneys were dried to constant mass and weighed to the nearest 0.001 g. Kidneys were finely ground in an acid-cleaned agate mortar. Concentrations of Cd, Cu, and Zn for the kidney samples (Apodemus sylvaticus, n = 157) were measured by ICP-AES with ultrasonic nebulisation for Cd, while Pb concentrations were measured by ICP-MS, both after total aqua regia digestion. Analytical quality was verified using blanks, CRMs (BCR 185 R, bovine liver; NIST 1547, peach leaves; DOLT-4, dogfish liver; DORM-3, fish protein) and duplicates (see Supplementary Tables S11 and S12 for details). . Before 2013, the capture of non-protected free-ranging rodent species immediately followed by their sacrifice was not considered to be an experiment on live animals and thus did not require protocol approval by an ethical committee. However, as stated above, care was taken to apply euthanasia protocols appropriate for small rodents in the field, and sampling authorizations were obtained from the appropriate administrative bodies.
Body condition and somatic indices. Body condition was assessed by the scaled mass index (SMI), and somatic indices were estimated using standard major axis (SMA) regression of ln-mass on ln-length as recommended by Peig and Green 25,26 . As the slope of this regression b SMA did not differ significantly between sites, it was estimated on the entire dataset excluding pregnant females for SMI (n = 154), and on the dataset excluding outliers for scaled somatic indices (n = 155). The SMA regression consists of estimating the predicted body mass (SMI) or the predicted organ mass (SLI for the liver, SKI for the kidneys) for each individual i when body length is standardised. Calculations for SMI follow the equation: where m i is the body mass and L i the body length of the individual i, b SMA is the slope of the regression of ln-body mass on ln-body length and L 0 the arithmetic mean of body length for the population. The scaled liver index SLI i and scaled kidney index SKI i values were similarly computed using organ mass instead of body mass. In this study, the b SMA values were 2.90 (95% confidence interval: 2.62-3.20) for SMI, 3.54 (95% confidence interval: 3.10-4.05) for SLI, and 2.91 (95% confidence interval: 2.55-3.31) for SKI. According to Green (2009, 2010), the b SMA value for SMI usually lies between 2.5 and 3.2, which can be used as a guideline to identify reliable estimates of the allometric exponent in mammals 28 .

Weighted TM concentrations in topsoils.
To account for the mobility of the wood mouse (home range of 2500 m²) 54 , a weighted TM concentration of the corresponding topsoil was calculated for each individual. This weighted concentration was calculated as the average of the TM concentration in the topsoil of the 8 plots surrounding the plot where the wood mouse was captured, with the capture plot being weighted twice. These weighted TM concentrations were used in the statistical models.
Data processing and statistical treatment. The Quantum GIS free software 55 was used for mapping.
Statistical treatment and regression analysis used the smart, lmodel2, stats, and pgirmess packages from the R software 56 . Topsoils. The non-parametric Kruskall-Wallis test was used to assess the differences in topsoil TM concentrations between sites in each park because the residues were not normally distributed. When the Kruskall-Wallis test was significant (p < 0.05), pairwise comparisons were made using Steel-Dwass post-hoc tests 57 .
TM in wood mice. A parametric ANOVA test was used to assess the differences in TM concentrations in wood mouse kidneys between sites in each park. When the ANOVA test was significant (p < 0.05), pairwise comparisons were made using Tukey's HSD test. Linear models were used to investigate the relationships between TM concentrations in kidneys and explanatory variables. Multivariate models were built with biological variables (mass and sex; age and size were dropped because of their relation to mass), geographical variable (site), and corresponding weighted TM concentrations in soils as explanatory variables. Biologically meaningful interactions (soil concentrations × sex and mass × sex) were also taken into account. As the range of TM concentration values in relation to sites did not always overlap, testing the interaction site × concentration in soil was not allowed 58 . Statistical treatments for Cd were only performed with the M1, M2, C1 sites, as M0, C0 and M2 presented too many censored data for Cd in soils. The best-fit model was selected using a backward stepwise regression. The drop1{stats} function combined with an F test was used; this method corresponds to a type II ANOVA. For each step, the explanatory variable with the largest p-value was deleted until the final step, giving the best-fit model, where all p-values were below the α = 0.05 level (Supplementary Table S13 for model selection). ANOVA tests were performed on the selected models. Model normality was examined by looking at plots of the standardised residuals versus leverage. Model outputs were satisfactory.
Body condition and somatic indices. Backward stepwise regression was also applied to determine whether body condition and somatic indices varied with geography, individual characteristics, or levels of individual contamination. Multivariate models were built with biological variables (only sex, as we considered SMIs), a geographical variable (site), and the four TM concentrations in kidneys as explanatory variables. Biologically meaningful interactions (kidney TM concentrations × sex) were also taken into account (Supplementary Table S14 for model selection).
Fluctuating asymmetry. In this study, six bilateral morphometric traits were selected for fluctuating asymmetry (FA) assessment: length and width of the three lower molars (Supplementary Fig. S5). All measurements were performed twice to control for measurement error. The FA consists of subtle random variations between each side (right, R and left, L) of bilateral traits that are supposed to be perfectly identical. These variations reflect the inability of individuals to correct errors occurring during early development. It has been shown that both genetic and environmental stresses decrease developmental stability 59 . The FA has been proposed as a useful tool to assess individual quality 30 . In fact, three types of biological asymmetry can be distinguished on the basis of the analysis of right minus left (R − L) frequency distribution: directional asymmetry (DA), antisymmetry (AS) and FA. The DA shows a pattern of normal R − L variation distributed about a mean point that is significantly different from zero. The AS shows a pattern of R − L variation distributed about a mean point of zero, but with a frequency distribution departing from normality 30 . The FA shows a normal distribution of R − L values with a mean of zero. Among these three asymmetries, only FA provides an estimation of developmental instability. The presence of DA and AS, together with measurement error, can bias FA estimation. A series of preliminary analyses was performed for each study trait as recommended by Palmer 29 . Individual FA levels were then estimated for each trait using absolute asymmetry. Linear models were computed to assess the relationship between |R − L| values and kidney TM concentrations. Population FA levels were estimated for each trait using FA10, i.e., between-sides variance corrected for measurement error, obtained from the results of linear mixed models with sides (fixed) × individuals (random) (see Palmer 29 , for details). Fisher tests were then performed for each trait studied to explore inter-site differences.
Data availability. All data generated or analysed during the current study are included in this published article and the Supplementary Information file.