The sampling sites are located near the Horne Smelter (Rouyn-Noranda, region of Abitibi-Témiscamingue, Québec, Canada) (Figure 1). Previous studies reported that the most contaminated sites were downwind of the prevailing winds (eastward elongation; Bonham-Carter et al. (2006)) and within the first 100 km from the smelter (Ettler, 2016; Telmer et al., 2004). Based on these results, sampling sites (n = 24) were selected along a 52 km transect starting from the Horne Smelter and oriented towards the southeast (Figure 1; Supplementary Table S1). All sites are far from other potential sources of TEs such as roads, landfills, settling ponds and mines. To minimize inter-site variability resulting from forest diversity near the Horne Smelter, sampling was only carried out in jack pine stands, an environment with rock outcrops where jack pine (Pinus banksiana), black spruce (Picea mariana), lowbush blueberry (Vaccinium augustifolium), and sheep laurel (Kalmia augustifolia) are commonly found (Kenkel, 1986). These conditions also meet the optimal ecological conditions for finding reindeer lichen C. rangiferina as well as wolf spiders Lycosidae (Koponen, 1987). In addition, the landscape topography of the specific sites is suitable for our purposes because they all have similar elevation (from 282 to 366 m) and they are directly exposed to the atmospheric flow carrying stack emissions loaded with TEs from the smelter (Aznar et al., 2008).

All lab wares (i.e., tubes, microtubes, spatulas) were soaked in 15% (v/v) nitric acid and rinsed with milli-Q water several times to minimize contamination. For PUF filters and lichens, a microwave-based digestion was performed (Gaga et al., 2019; Charette et al., 2021), while for spiders, a less aggressive digestion using an incubator was applied (Rosabal et al., 2012). Each PUF filter was divided into thirds of ∼100 mg each for determination of total TE concentrations. Microwave digestion of these thirds was performed in a mixture of nitric acid (7.5 mL; 67%–69% HNO₃, Optima grade, Fisher Scientific) and hydrogen peroxide (2.5 mL; 30% H₂O₂, Certified ACS, Fisher Scientific) for 20 min at 200°C in sealed polytetrafluoroethylene (PTFE) vessels (Multiwave 5000, Anton Paar) (Gaga et al., 2019). Ultrapure milli-Q water was then added to reach a total volume of 100 mL before measurements. C. rangiferina samples were freeze-dried for 24 h and homogenized (reduced to fine dust). Analysis of total mercury (THg) was performed on freeze-dried lichen samples using a direct mercury analyzer (DMA-80, Milestone Corp.). Digestion of lichen samples for determination of total TE concentrations was performed in a solution of hydrochloric acid (1 mL; 32%–35% HCl, Optima grade, Fisher Scientific), nitric acid (3 mL; 67%–69% HNO_3, Optima grade, Fisher Scientific) and hydrogen peroxide (500 μL, 30%–32% H₂O_2, Optima grade, Fisher Scientific). To complete the digestion, samples were heated to 200°C for 20 min in closed pressurized vessels in the microwave oven (Multiwave 5000, Anton Paar) (Charette et al., 2021). Ultrapure milli-Q water was then added to reach a total volume of 15 mL. Lycosidae samples were freeze-dried for 24 h and homogenized before digestion in nitric acid (500 μL, 67%–69% HNO_3, Optima grade, Fisher Scientific). They were heated at 65°C in an incubator for 4h. Once samples cooled, hydrogen peroxide (250 μL; 30%–32% H₂O₂) was added to oxidize organic matter (Rosabal et al., 2012). Ultrapure milli-Q water was then added to reach a total volume of 15 mL before TE measurement. All TE measurements were performed using an inductively coupled plasma-tandem mass spectrometer (ICP-QQQ, 8900 Triple Quadrupole, Agilent). Method detection limit for each studied TEs is shown in Supplementary Table S3.

To ensure the accuracy of measured concentrations and to assess the degree of TE extraction of our different digestion protocols, various certified reference materials were also analyzed. Reference materials NIST1573a (tomato leaves, National Institute of Standards and Technology, USA) and BCR-670 (Duck weed, Institute for Reference Materials and Measurements, European Commission) were used to assure quality measurement in lichens for conventional TEs and REEs, respectively. For Hg measurements, blanks were analyzed every 10 samples (No Hg contamination was observed), and certified reference materials (CRMs) TORT-3 (lobster hepatopancreas, National Research Council of Canada) and DORM-2 (fish protein, National Research Council of Canada) were used for quality control. Average Hg recoveries were 97.8% ± 0.3% (n = 3) for DORM-2 and 97.8% ± 2.4% (n = 9) for TORT-3. For spiders, the digestion was evaluated using DOLT-5 (dogfish liver, National Research Council of Canada) for conventional TEs and BCR-668 (mussel tissue, Institute for Reference Materials and Measurements, European Commission) for REEs as these CRMs have already been used for invertebrates (Marginson et al., 2023). The recovery percentages (Supplementary Table S3) for most TEs under study varied from 70% to 120%, with the notable exception of La in BCR-668 which did not meet our quality control requirements; La data were therefore removed from this data set. The detection limit for each TE measured is shown in Supplementary Table S3. NIST1573a, BCR-670 and BCR-668 were used because they are the only CRMs available for REEs in the analyzed matrices.

All numerical data are presented in dry weight and as the arithmetic mean (± standard deviation, SD). Statistical tests were performed using SigmaPlot (10.0) and R (3.6.1). Normality was confirmed using Q-Q plots and the Shapiro-Wilk test. Outliers were identified using Grubbs test and removed from the dataset (<1% of samples). Regression analyses were performed to explore the relationship between TE concentrations and distance from the Horne Smelter. Exponential decay relationship was explored for all metals and three parameters (y ₀ , a, b) were determined (Equation 1) when significant regressions were obtained. The determination coefficient (R²), regression line and equation were reported when the model was significant (p < 0.05). The y ₀ parameter of Equation 1 was used to estimate TE baseline concentrations. y = y 0 + a · e − b x (1)

Spatial metal concentration gradients ([M]_maximum/[M]_minimum) for each monitor collected over 52 km were calculated using mean TE concentrations in each sampling site. To determine spatial differences compared to the site closest to the Horne Smelter (i.e., the presumed most contaminated site), a Kruskal-Wallis test followed by Dunn’s test (after Bonferroni correction) was applied for TE measurements in samples of PAS, lichen and spiders. The sum of REEs (ΣREE) included Y (and excluded La) as it was strongly correlated to the other REEs (Supplementary Table S4; R² = 0.85–0.99). REE concentrations in lichen were normalized by the Post-Archean Australian Shale values (PAAS) to assess possible anomalies in their relative proportion, which were calculated using Equation 2 (example from Ce) (Dang et al., 2023; Pourmand et al., 2012). C e a n o m a l y = 2 · C e P A A S L a P A A S + Pr P A A S (2)

Non-linear models were used to describe the relationship between C. rangiferina ground cover (%) and distance from the smelter Equation 3. Also, the exponential decay relationship (Equation 1) was explored to relate TE concentrations in lichens with the C. rangiferina ground cover (%). y = a 1 + e − x − b c (3)

The concentrations in deployed PUF-PAS showed significant exponential decay relationships with the distance from the Horne Smelter (p < 0.05) for a variety of TEs, including As, Ba, Cd, Cu, Fe, Mn, Mo, Pb, V, Zn and some REEs (Y, Pr and Nd). The negative relationship for Co was linear (p = 0.003) (Figure 2; Supplementary Figure S1). The exponential decay relationships observed between TE concentrations in filters and the distance from the smelter showed a high level of association (based on R² values), varying from 0.50 (Zn, Y) to 0.97 (Pb, Se) (Figure 2; Supplementary Figure S1). High spatial concentration gradients (ratio: [M]_maximum/[M]_minimum) were observed along the 52 km transect sampled, where this ratio ranged from 1.7 (Co), 4.5 (Zn), 83.3 (As) to 375 (Cu) (Supplementary Table S5). For most metals, the PUF-PAS deployed after 28.7 km showed significant differences in TE concentration compared with the measurement in the PUF-PAS closest to the Horne Smelter at 1.4 km (i.e., the most contaminated site) (Supplementary Table S6).

Similar to the passive air sampler results, highly significant exponential decay relationships were observed between TE concentrations in biomonitor samples and distance from the Horne Smelter (Figures 3, 4, Supplementary Figures S2–S4). However, the decay relationship was much less steep compared to PUF-PAS. For lichens, this relationship was significant for a wide range of TEs (As, Ag, Ba, Pb, Se, Sb, Cu, Ti, Co, Cd, Hg, Ni, Zn, Fe, Mo, V, La, Ce, Pr, Nd, Sm, Gd, Y, Tb, Dy, Ho, Er, Yb) (Figures 3, 4;Supplementary Figures S1–S4), while for spiders, significant exponential decay relationships were observed only for As, Se, Ag and Pb (Figure 3). For Cu and Cd in spiders, a linear relationship with the distance from the Horne Smelter was observed. The spider samples showed TE concentrations that are, in general, considerably greater than the concentrations measured in lichens (more than 50 times higher for Cd). Based on concentrations in lichens and spiders, TE enrichment generally extended up to 28.7 km from the smelter (similar to PAS measurements), after which the bioaccumulation concentrations significantly differed from the TE concentrations of the more contaminated site located at 1.4 km (Supplementary Tables S7–S9). We observed high spatial bioaccumulation gradients for elements previously reported in the area (e.g., As, Cd, Cu, Zn, Pb) as well as for unstudied elements (e.g., Co, Mo, V, Ag, REEs) (Supplementary Table S10).

Despite clear REE enrichment in lichens, where Nd showed the strongest spatial gradient among REE members (19.3), the REE patterns in lichen samples reflect their relative natural abundance in bedrock (Supplementary Figures S5, S6). The average REE concentrations in lichens display a log-linear pattern decreasing with atomic numbers, following the Oddo-Harkins rule, and is consistent regardless of REE concentration (Supplementary Figure S6A). The PAAS-normalized REE concentrations in lichens showed that heavy REEs (i.e., Tb, Dy, Ho, Er, Yb) are slightly enriched compared to light REEs (i.e., Ce, Pr, Nd, Sm) (Supplementary Figure S6B). The REE patterns in lichen samples also exhibited a slight positive Gd anomaly that is consistent regardless of REE concentration.

We recorded the presence of C. rangiferina at a minimum distance of 2.7 km from the smelter. The site closest to the smelter was also surveyed but C. rangiferina was absent. Expressing the impact of reported TE contamination on C. rangiferina abundance, we observed that vegetated ground cover increased with the distance from the Horne Smelter. In the first 5 km from the smelter, C. rangiferina presence was sparse, with ground cover <1%. The abundance of lichens is considerably greater beyond 10–20 km from the smelter, although there was important variability in lichen ground cover within sites (Figure 5A). The bioaccumulated concentrations of Hg (p: 0.05; R² = 0.69) and Ag (p: 0.05; R² = 0.58; data not shown) were correlated with lichen abundance in the field (Figure 5B).

Airborne particulates emitted from smelting operations are often associated with higher concentrations of TEs such as As, Cd, Cu and Pb in the surrounding environments than the natural background values (Kasongo et al., 2024). Smelter-derived particulate emissions from a point source have been shown to induce halo-like contamination patterns, resulting from differential settling due to particle size (Csavina et al., 2011; Csavina et al., 2014; Berryman et al., 2024). The PUF-PAS measurements in this study are consistent with emissions from a point source and showed a clear enrichment within the first 28.7 km of the smelter for most of the TEs studied, which could be related to airborne particles emitted from the smelter. The highest maximum/minimum concentration ratios estimated were for Cu (375), Pb (184), and As (83), respectively, which is comparable with reported stack emissions and field observations in the vicinity of other copper smelters (Fernández-Camacho et al., 2012; Svoboda et al., 2000; Zhang et al., 2022). In addition to Cu, Pb and As, there were also significant spatial concentration gradients for TEs previously unreported in this region such as Ba, Mo, Y, Nd and Pr. The ability of PUF-PAS to detect REEs (i.e., Y, Nd, and Pr) in ambient air is supported by previous studies where dispersion of REEs near smelting facilities was found to be primarily mediated by the atmospheric deposition of fine REE-laden particles (Brewer et al., 2022; Wang et al., 2014). These newly detected elements might be associated with the smelter’s recycling of e-wastes and hazardous materials, suggesting these recycling actions contribute to environmental TE contamination (Bi et al., 2010; Song and Li, 2014). As such activities are expected to increase in coming years, more environmental monitoring efforts should be performed in the future for these emerging TEs.

According to our passive samplers, sites located further away in forested ecosystems (more than 28.7 km from the smelter) showed a stable and relatively low concentration of TE-rich airborne particles, which contrasted with the deposition of TE-enriched smelter dust found up to 100 km away from the Horne Smelter in previous studies (Knight and Henderson, 2006; Telmer et al., 2006). The implementation of emission reduction measures taken in recent decades may partly explain why PUF-PAS measurements detected atmospheric TE enrichment at substantially shorter distances compared to previous studies (Crouse et al., 2019; Diener and Mudu, 2021; Wolf et al., 2020). However, the magnitude of enrichment measured by the PUF-PAS should be viewed with caution as they showed variations in their affinity for different TEs and particle sizes (Gaga et al., 2019; Li et al., 2018). Metals that tend to be found in very fine particles (Cu, Ni) were shown to have significantly higher absorption rates in the PUF than those that are bound to the coarse fraction. In addition, it has been shown that the “unfilterable” phase, comprising gaseous species and solutes in filter-penetrated droplets, accounts for a considerable proportion of the atmospheric TE content, leading to underestimation when considering only the particle collectors (Zhang et al., 2022). Despite the fact that smelter emissions are notoriously subject to batch processes (Ministère de l'Environnement et lutte contre les Changements Climatiques, 2018), the PUF-PAS devices used captured a 3-month snapshot of TE contamination in airborne particles. Passive samplers have the potential to provide low-cost alternatives to traditional air sampling techniques and here they demonstrated their potential to meet the broad sampling requirements posed by smelters, particularly in remote or logistically constrained environments (Kasongo, et al., 2024).

Similar to PUF-PAS measurements, regression analyses of TE concentrations in both biomonitoring organisms showed contamination patterns typical of a point source emission (Ettler, 2016; Fry et al., 2020; Stafilov et al., 2010). Numerous studies have reported analogous patterns using terrestrial biomonitors such as lichens (Bari et al., 2001; Cloquet et al., 2006; Fry et al., 2020; Grodzińska et al., 1999; Gačnik et al., 2024), but none have measured a range of TEs (including contaminants of emerging concern) as broad as in this study. We found a high TE concentration gradient for Pb, Cd, Ni, Sb, Zn, and Cu, which is comparable with previous reports of particulate contamination and likely reflects the major TEs emitted from the Horne Smelter. Consistent with the atmospheric particles analyzed from the PAS, various elements (e.g., Mo, REEs) were newly reported for this region, showing relatively slight TE gradients, but with significant spatial differences compared to the TE concentrations of the closest site to the Horne Smelter. These results reveal gaps in the monitoring of smelter-derived inorganic pollution in the Rouyn-Noranda region, particularly for emerging contaminants such as REEs, which are currently unregulated (Gwenzi et al., 2018). Since several TE measured in lichens are considered strategic and critical metals (e.g., Cu, Ni, Co, REEs, etc.) because of their increasing use in electronics, a link between their emission and the recycling of e-waste and hazardous materials at the smelter is suspected (Song and Li, 2014). Likewise, Sb, also found in the lichens in this study, has been shown to be directly related to emissions from e-waste recycling facilities (Bi et al., 2010). The presence of V in lichens suggests that some particles may be emitted from the sulfuric-acid plant adjacent to the smelter, where V₂O₅ is used as a catalyst to oxidize SO₂ to SO₃ (Samuel et al., 1989; Zdanowicz et al., 2006).

Lichen samples also showed TE enrichment extended to approximately 28.7 km from the smelter in the southeast direction, but still lower than the spatial TE deposition reported mostly 20 years before (Knight and Henderson, 2006; Telmer et al., 2006). Although PUF-PAS captures fine atmospheric particles and lichens are able to trap coarse particles within their medulla, both monitors were jointly able to detect significant metal concentration differences in the last sampled 35 km compared to the closest sites (i.e., 3-35 km) to the Horne Smelter. This consistency could reflect a high connection between atmospheric emissions and deposition of various TEs emitted from the Horne Smelter in sampled forest ecosystems. Interestingly, these spatial differences are not only observed for TEs previously found in the area (i.e., As, Cd, Cu, Pb) but also for newly detected TEs (REEs, Co, Mo). In addition, the fact that these two monitors, which provide different time-integrated information (for example, 3 months for PUF-PAS, some years for lichen), yielded similar transect spatial results, means no significant changes in TE deposition patterns occurred in recent years.

Regarding the magnitude of TE accumulation in the vicinity of a smelter, Grodzińska et al. (1999) used the lichen Cladonia stellaris to conduct a biomonitoring study near the Ni-Cu smelter complex in Monchegorsk, one of the most polluted cities in Russia. Compared to their results, in our study lichens sampled at the site closest to the Horne Smelter had similar TE concentrations for Fe, Al, Cu, Cd, Cr, Pb and Zn, but an order of magnitude lower for Co, Mn and Ni. Similarly, a study conducted near the Cu-Ni smelter in the Finnish city Harjavalta reported concentrations of Cu, Ni, Zn, Fe, Mn, Cd, and Pb in C. rangiferina that aligned closely with the maximum concentrations observed in our lichens (Salemaa et al., 2004) (Supplementary Table S11).

Metal concentrations for all studied TEs (except Cr) in areas farther from the smelter were estimated by exploiting the relationship between their concentration in lichens and the distance from an isolated point emission source (Supplementary Table S12). A large-scale study of the TE concentrations in peltigeralean lichens in northern Québec showed Al, Ti, V, Fe, Co, Ni, Cu, Zn, Mo, and Cd levels similar to those calculated from our regression-based estimates of the baseline values for C. rangiferina lichens (Darnajoux et al., 2015). However, our regression-based approach showed Mn concentrations that were an order of magnitude lower, and Pb levels that were 6 times higher than those from Darnajoux et al. (2015). Furthermore, estimated C. rangiferina Cd, Cu and Pb concentrations exceeded those found in terricolous lichen Peltigera canina from Yukon Territory, Canada (Pouillé et al. 2024). It is worth noting that comparisons between species of differing genera may be questionable due to the specific physiology and distinct morphologies of the species. Fortunately, a study conducted by Chiarenzelli et al. (2001) in a remote region of northern Canada (near Otter Lake, Nunavut) showed that concentrations obtained using C. rangiferina are of the same level as our measurements for most TEs, including As, Ag and REEs.

While metal quantification in C. rangiferina provided valuable data, it remains unclear to what extent the observed TE enrichment is due to current smelter operations and to what extent it is due to resuspension of soil particles containing historical TE deposits. Widory et al. (2018) used the relationship between ²⁰⁶Pb/²⁰⁴Pb isotope ratio and Pb concentration in C. rangiferina samples to determine its main atmospheric sources in the Rouyn-Noranda area. Their results showed that Pb enrichment east of the city is the result of binary mixing between emissions from the Horne Smelter and regional Pb background. Although a rigorous source apportionment analysis was not the focus of this study, the Pb enrichment in lichens showing a strong correlation (R² ≥ 0.89) with the concentrations of Ni, V, Mo, Ag, Sb, Cu, As, Co, Se, Cd, Fe, Al, Zn, ΣREE across 2 orders of magnitude suggests a similarly proportioned mixed origin for these TEs (Supplementary Table S13). A slightly weaker correlation for Mn (R² = 0.83) compared to the other TEs could be due to a relatively higher proportion of Mn from geogenic sources (Bargagli et al., 2002). As for Ba (R² = 0.75), Ti (R² = 0.63) and Cr (R² = 0.51), the weaker association with Pb concentrations in lichens may be related to a more pronounced mixing with traffic-related emissions (Birmili et al., 2006; Sakata et al., 2021).

The REE log-linear decreasing pattern observed in lichen samples is similar to that previously observed in other sites (MacMillan et al., 2017; Picone et al., 2022) (Supplementary Figure S6A). However, PAAS-normalized concentrations in lichens typically show a clear downward slope, with higher light REEs enrichment relative to heavy REEs. Our results instead showed a somewhat small deviation from a horizontal line, with even a slight enrichment in heavy REEs (Supplementary Figure S6B). This trend should be interpreted with caution because the concentrations further from the suspected source were very low. Apart from the slight Gd anomaly, which may be associated with anthropogenic sources (Lafrenière et al., 2023), our results essentially indicate that the bioaccumulation patterns were relatively uniform across the REE group. This uniformity leads us to suggest that REE enrichment in the vicinity of the Horne Smelter is not particularly associated with e-waste and hazardous recycling activities, but rather to complex concentrate inputs to the smelter containing some amount of REEs. However, analyzes of representative samples of the smelter feed materials would be required to conclusively establish this, especially since no REE enrichment has ever been identified near a copper smelter with recycling activities (Parviainen et al., 2019).

Measurements of C. rangiferina ground cover shows that this biomonitor does not grow near the Horne Smelter. The same situation is observed near the Harjavalta smelter, where C. rangiferina is only present sparsely starting from 2 km away (Salemaa et al., 2004). The sigmoidal relationship between distance from the smelter and abundance of C. rangiferina at the sampling sites is consistent with what would be expected near a point contamination source. Significant correlations between the abundance of C. rangiferina in the field and its TE bioaccumulation were observed, particularly for Hg (R² = 0.69) and, although to a lesser extent, for Ag (R² = 0.58). These findings are consistent with the results of Puckett (1976) who reported that both elements are the most toxic metals to lichens, altering their cellular membrane and inhibiting photosynthesis.

With respect to spiders as biomonitors, spatial differences were observed only for As, Se, Cu, and Cd. Correlations have been shown between TE bioaccumulation in invertebrates and the soil TE content (Heikens et al., 2001), indicating the availability of these elements from smelter-contaminated soils. Mean concentrations of bioaccumulated Cu and Cd in spiders from the sites closer to the Vaudray-Joannès Lakes Biodiversity Reserve are at comparable levels, and of the same order of magnitude as studies in other part of the world using ground-dwelling spiders as biomonitors (Jung et al., 2005; Larsen et al., 1994). Previous research conducted near the Vaudray-Joannès Lakes Biodiversity Reserve showed limited inter-individual variability of TE concentrations in Lycosidae (Ponton et al., 2018), but our measurements were more variable within sites, limiting the observation of spatial trends along the studied transect for all TEs. The spider samples for this study were only standardized for weight, which may explain the amount of interindividual variation observed. It is well known that biotic factors such as diet, sex, and age structure can affect the bioaccumulation of TEs in organisms. Also, the lack of spatial differences for the other TEs studied compared to lichens (which is a primary producer) indicates that spider exposure to polymetallic contamination is likely associated, in part, with the dietary pathway as a consumer, rather than being solely atmospheric like lichens. Spiders are also able to reflect the TE deposition from the smelter, but some ecological processes may alter this spatial TE concentration patterns.