Introduction

Among anthropogenic activities, metal mining and metallurgic industry have led to extreme environmental changes through pollution of large areas with Trace Metal Elements (TME) like zinc (Zn), cadmium (Cd) or lead (Pb). Soil TME concentrations in so-called metalliferous environments are often so toxic that only few living organisms are able to survive and reproduce. This is particularly true for plant species that are sessile organisms that cannot avoid toxicity. Some of them, called metallophytes, nevertheless evolved specific mechanisms of TME tolerance that allow colonizing metalliferous environments (Antonovics et al. 1971).

Some TME tolerant species also evolved unusual abilities of TME accumulation in aboveground parts, without suffering toxicity symptoms. Above an accumulation threshold that differs among metals and takes into account the environmental phenotypic variance, these plants are called hyperaccumulators (van der Ent et al. 2013). Metal hyperaccumulation is known as one of the most intriguing traits in evolutionary biology. However, its primary adaptive function, as well as the current selective factors that make it evolve, are still debated (Pollard et al. 2014; Rascio and Navari-Izzo 2011). Several non-mutually exclusive evolutionary hypotheses were previously suggested (Boyd and Martens 1992, 1994, 1998). Hyperaccumulation could have served as a strategy of TME tolerance, elemental defense against herbivores or pathogens, elemental allelopathy against competitors, drought resistance, or may be an inadvertent consequence of some nutritional processes. In addition, three hypotheses, also non-mutually exclusive, have been specifically suggested to explain the existence of “facultative hyperaccumulators”, i.e. hyperaccumulators that are facultative metallophytes (syn. Pseudometallophytes) that display species-wide but not necessarily uniform hyperaccumulation capacities. The “phylogenetic hypothesis” states that ancestors evolved on metalliferous environments, and that, after expansion of species range, non-metallicolous individuals conserved hyperaccumulation capacities. The “incremental advantage” hypothesis postulates that hyperaccumulation originally conferred some adaptive advantages on non-metalliferous soil regarding multiple biotic and abiotic selective factors, and thus can be maladaptive on metalliferous soil unless specific detoxification mechanisms secondarily evolve. The “inadvertent uptake” hypothesis suggests, as previously raised by Boyd and Martens (1992) for hyperaccumulators in general, that hyperaccumulation may be a side effect of some other nutritional processes. In the latter hypothesis, hyperaccumulation would then be non-adaptive or even maladaptive.

In summary, competing hypotheses consider that hyperaccumulation could have evolved either in metalliferous and non-metalliferous environments, and could be either neutral or adaptive, evolving as a response to either abiotic or biotic selective pressures. However, with the exception of the elemental defense hypothesis that benefits from an abundant literature (Hörger et al. 2013), the debate usually remains mostly theoretical because alternative hypotheses have never been definitely refuted (Boyd 2007). This may be due to the fact that “abiotic” hypotheses do not clearly state the nature of selective pressures that could have acted on hyperaccumulation.

Identifying the selective pressures that influence the evolution of an adaptive trait is a challenging issue in evolutionary biology. When there are strong presumptions about the role of one biotic or abiotic parameter, its effect may be directly estimated in controlled conditions, using for example experimental evolution approaches. Otherwise, indirect approaches consist in testing, in natural conditions, the correlation between phenotypic variation of the trait under selection and environmental variation of any parameter possibly representing a selective pressure (Linhart and Grant 1996; Merilä et al. 2001). This approach, based on statistical correlations, is not a conclusive demonstration. However, it might provide some insights into putative selective pressures, whose effect can be observed on several generations of selection in natural populations.

In search for the “raison d’être” of hyperaccumulation, facultative hyperaccumulator species may be the more appropriate models because they present both metallicolous and non-metallicolous populations and allow following metal hyperaccumulation in various conditions of metal exposure. Arabidopsis halleri (Brassicaceae) is a model pseudometallophyte that tolerates and hyperaccumulates Zn and Cd. In controlled conditions, A. halleri displays species-wide Zn and Cd hyperaccumulation capacities, despite among-population quantitative variation of hyperaccumulation levels (Bert et al. 2000, 2002; Macnair 2002). Such variation is a prerequisite for adaptive evolution. Besides, it has been demonstrated that its population structure was not random and may result from non-neutral evolution. Indeed, shoot Zn concentrations estimated in experimental conditions were on average higher in non-metallicolous populations than in metallicolous ones (Bert et al. 2002). In field conditions, however, shoot Zn concentrations are expected to be lower in non-polluted conditions than in polluted conditions, due to the lower levels of available zinc in soils. This discrepancy between observations in field and experimental conditions confirms that the determinism of hyperaccumulation levels is not fully genetic and must involve interactions between genotype and environment (as demonstrated for Noccaea caerulescens, the other model pseudometallophyte: Escarré et al. 2013). It also highlights that difference may exist among individuals in their capacity to translocate metals in shoots. Indeed, the plant/soil concentration ratio, also called the bioconcentration factor (Lovy et al. 2013), is usually expected to be on average higher in non-metallicolous individuals, as demonstrated in situ for N. caerulescens (Molitor et al. 2005) and also recently for A. halleri (Stein et al. 2017).

In this context, this study aimed at investigating the potential role of environmental abiotic parameters in driving the microevolution of hyperaccumulation in a set of neighboring sampling sites of A. halleri from the Italian Alps. To achieve this objective, we characterized various habitats where A. halleri occurs, in a landscape expected to show environmental heterogeneity in particular for metal concentrations in soils. Indeed, several spots of metallurgic activities are scattered in this region of Italy. In this study, A. halleri habitats were described as extensively as possible in terms of abiotic parameters (topography, climate, elemental soil concentrations). Floristic data were also collected, as plant communities give integrative information about the environmental conditions. Finally, we tested for correlations between spatial variation of abiotic parameters and natural variation of the Zn and Cd accumulation levels. For this purpose, hyperaccumulation abilities were estimated both in terms of shoot metal concentrations and bioconcentration factors (hereafter referred as “hyperaccumulation variables”). Plant metal concentrations were assessed in field conditions, and soil metal concentrations were directly measured in the rooting soil of the target plants in order to estimate field bioconcentration factors.

Materials and methods

Description of sampling sites

In spring 2008 and 2009, two valleys in Lombardy (Italy) have been explored (Fig. 1). 17 locations have been sampled in the Val del Riso valley, where A. halleri occurs in diverse habitats (Table 1). The Val del Riso valley notably hosted mining and smelting activities at high and low altitude, respectively, between 1850’s and 1980’s. From these metallurgic sites, soil deposits and atmospheric fallouts can be expected to have caused gradients of metal concentrations in soils. In the Val del Riso valley, we thus distinguished 9 so-called “M” sites that were located in mining or smelting areas from 8 so-called “NM1” sites that were outside these areas. In addition, 6 locations have also been sampled in a second valley 40 km away from the first one, where no mining or industrial activity was reported (Fig. 1). For this reason, these distant sites were called “NM2”. NM2 sites included two sites on clay soils along a road (Val Paisco) and four sites on calcareous soils near a small village (Sommaprada) (Table 1).

Fig. 1
figure 1

Geographic localisation of the Italian Arabidopsis halleri sampling sites. Populations I12 to I22, I24 to I27 and I35-I36 are located in the “Val del Riso” valley, populations I28-I29 are located in the “Val Paisco” valley, and populations I30 to I34 are located in the “Sommaprada” valley. Red stars: metalliferous sites; black points: non-metalliferous sites

Full size image
Table 1 General characteristics of Arabidopsis halleri Italian sampling sites
Full size table

Field work: Topographic, climatic and floristic data

Altitude, slope, azimuth and absolute humidity of air were directly recorded in the 23 sampling sites, providing one single value by sampling site. Absolute humidity, expressed in gram of water by kilogram of air, was measured with a thermo-hygrometer HD200 (KIMO Instruments®). Because it is based on mass of air, this measure is not sensitive to changes of air volume that can occur under various temperatures or pressures. Then, for each sampling site, the maximum daily direct solar irradiance (kWh m−2 d−1) was calculated from latitude, slope, and azimuth data. Note that this calculation only gives the theoretical direct solar irradiance reaching a given surface, since it assumes that there is no shading effect by surrounding objects (like buildings or trees) and a clear-sky situation. In addition, in order to characterize A. halleri habitats in terms of vegetation and successional stages, the occurrence of the most abundant accompanying herbaceous species was recorded. We checked plant species taxonomy using TAXREF (Gargominy et al. 2016). We described habitat types following the European Nature Information System classification (EUNIS, http://eunis.eea.europa.eu/), on the basis of dissimilarities among sites, performed on plant communities (see Statistical analyses).

Soil and plant analyses

In 21 of the 23 sampling sites (I20 and I26 not included), 5 rosettes and their corresponding rooting (30 cm depth) soil samples were collected every 5 m. Plant samples were cleaned to remove any soil particle and then air-dried at 55 °C during 72 h. Soil samples were air-dried at 40 °C during at least 3 weeks and sieved with a 2 mm mesh for pH and elemental analyses. Soil pH was measured with a glass electrode on a water-saturated sample (2.5 g of sieved soil in 10 mL of distilled water) and after stirring for 20 min. Soil mineral elements were extracted with ammonium acetate-EDTA 1 N (pH 4.65) for 30 min (10 g soil in 50 mL). Element concentrations of Ca, Cd, Cu, Fe, K, Mg, Na, P, Pb and Zn were determined by Induction Coupled Plasma-Optical Emission Spectrometer (ICP-OES, Varian Vista MPX). Dried and mixed plant samples were mineralised in a mixture of nitric and perchloric acid (2:1) and then analysed by ICP-OES for the same elements as in soil. However, only Zn and Cd shoot concentrations were included in statistical analyses because A. halleri is Zn and Cd hyperaccumulator. In order to respect technical limits of ICP-OES and thus retain only reliable values, plant samples whose weight was less than 50 mg were discarded from the final dataset.

In 22 of the 23 sampling sites (i.e. all but I21), soil particle size distribution was estimated by successive sievings (INRA of Arras, France) to determine the rates of stones (>0.5 cm), gravels (0.2 to 0.5 cm) and fine particles (<2 mm) from 3 samples among the 5 collected in the field.

Statistical analyses

In order to explore the environmental data and examine correlations between abiotic parameters, we performed a PCA analysis using the Factoshiny package built for R 3.3.3 (Vaissie et al. 2015). We also performed a Non-Metric Multidimensional Scaling (NMDS) analysis on floristic data to complete the PCA analysis based on abiotic parameters and thus fully characterize A. halleri habitats (package Vegan, R 3.3.3). This analysis provides a 2-D visualisation of the similarities among sites based on the species composition of the communities. NMDS analysis was based on the Jaccard’s similarity index. Jaccard’s index is designed for the analysis of species occurrence data and meets the asymmetry criterion usually preferred in ecological comparisons of communities, i.e. we abstain from drawing any ecological conclusion from the absence of a species in two sites (Legendre and Legendre 1998).

The environmental parameters as well as the four hyperaccumulation variables (Zn and Cd shoot concentrations and bioconcentration factors) were also individually treated to test for differences between metallurgic and non-metallurgic areas (M vs NM1/2) and between valleys (NM1 vs NM2). For this purpose, non-parametric Kruskal-Wallis tests for 3 independent groups (M, NM1 and NM2) were performed, followed by pairwise comparisons using Tukey and Kramer (Nemenyi) test. In the latter analyses, depending on the number of values recorded in each sampling sites, sample sizes were [M = 32; NM1 = 21; NM2 = 19] for hyperaccumulation variables and soil element concentrations, and [M = 9; NM1 = 8; NM2 = 6] for other abiotic factors. All these analyses were implemented with R software (R 3.3.3).

To investigate to what extent some abiotic parameters could explain variations of Zn and Cd shoot concentrations and bioconcentration factors, Partial Least Square Regressions (PLSR) were implemented for each hyperaccumulation variable. PLSR provides combinations of explanatory variables that maximize the variation explained in a biological response variable such as seed dormancy (Wagmann et al. 2012) or phenological traits (Brachi et al. 2013). PLSR is particularly useful when the number of observations is inferior to the number of explanatory variables, and when there is multicollinearity among explanatory variables (Tenenhaus 1998). PLSR were conducted for all data and then separately for M, NM1 and NM2 groups, leading to a total of 16 models. Sampling sites that could not be analysed for several abiotic factors (I20, I21 and I26) were excluded from the regression analysis. When the number of soil samples was lower than the number of plant samples (for plant growing on a too thin layer of soil), missing data were replaced by mean values of the sampling sites. Finally, because only one single value per population was available for topographic and climatic data, each value was replicated for all samples from a same sampling site. In the final dataset, between 1 and 5 pairs of plant/soil samples per sampling site were available for PLS analyses, with a mean of 3.4 pairs. Regression analyses were implemented with the PLS procedure of SAS (SAS Institute Inc 2011), using the “simpls” method. For each analysis, the number of extracted factors was set by “leave-one-out” cross validation and the minimum value of the root mean Predicted Residual Sum of Squares (PRESS). For each regression model, the most significant environmental parameters displayed “Variable Importance in Projection” (VIP) values higher than 0.8. Effect direction and intensity of each parameter were respectively specified by the sign and the absolute value of the corresponding centered and scaled regression coefficients.

Results

Raw data on plants collected in the field, corresponding soil analysis, and topographic and climatic measurements are presented for the 23 sampling sites of Arabidopsis halleri in Table S1. Presence data of the different taxa recorded in the same sampling sites are presented in Table S2.

General habitat characterization using abiotic data

The two first axes of PCA both explained 47.56%, with 32.39% for the first axis and 15.17% for the second one (Fig. 2). Overall, M sampling sites were largely distributed over the graph of individuals, whereas NM1 and NM2 sampling sites aggregated in different areas of the graph (Fig. 2a).The first axis clearly separated NM sites from some M sites (I16, I19, I12, I15 and I35, contrary to I13, I17, I18 and I35: Fig. 2a and b) mainly because of higher Cd, Ca, Cu, Pb and Zn soil concentrations, lower Na soil concentrations, and coarser soil structure in some M soils (Fig. 2c). According to the first axis, the I13, I17, I18 and I36 M sites would be somehow comparable to NM sites regarding their abiotic parameters (Fig. 2b). However, we are aware that sampling sites were represented by too few points to be fully characterized, and that our results mainly reflect general differences among M, NM1 and NM2 groups. The second axis tended to split NM sites into NM1 and NM2 groups (Fig. 2a) mainly because of higher Mg soil concentrations and absolute humidity values in NM1 sites, and higher Fe soil concentrations, solar irradiance values and altitudes in NM2 sites (Fig. 2a and c).

Fig. 2
figure 2

PCA results related to the abiotic parameters recorded in the field, showing (a) the graph of the individuals identified according to their group (M, in black; NM1, in red; NM1, in green), or (b) to their population of origin (in black colour). The graph of variables (c) gives the correlations among the abiotic parameters. The elemental soil concentrations are indicated using the abbreviation of the element (Ca, Cu, Cd, Fe, K, Mg, Na, P, Pb, Zn). irr: maximum solar irradiance; hum: absolute humidity; stone: rate of stone; grav: rate of gravels; fine: rate of fine particles; alt: altitude; pH: soil pH

Full size image

General habitat characterization using floristic data

Prior to perform analyses, Taraxacum species have been grouped into one single species complex. Dubious species were named as Genera sp. when found into one site only, but evicted when found in more than one site, due to the risk of merging several species from different habitats.

The stress value of the NMDS (the disagreement between 2-D configuration and predicted values from the regression) was 0.146, and the stressplot (the obtained ordination distance against the observed dissimilarity) was almost increasing monotonically (linear fit R2 = 0.91, non-metric fit R2 = 0.98). These values suggest that the NMDS provided a great representation of the data in reduced dimensions.

While M sites were relatively spread over the NMDS graph, NM sites (Fig. 3, in red and green) were grouped altogether, suggesting that NM sites are more similar in terms of species composition than M sites. Among NM sites, I14, I20, I22, I26 and I30–32 were densely aggregated. I14, I20, I22 and I26 (NM1) were mainly associated with plant species occurring in oligo to mesic calcareous meadows in their low-altitude to subalpine forms (Astrantia major, Anthoxanthum odoratum, Carum carvi, Galium mollugo, Geranium sylvaticum, G. phaeum, Heracleum sphondylium, Plantago lanceolata, Trifolium pratense, Phyteuma orbiculare, Ranunculus acris, R. repens, Rumex acetosella – EUNIS habitat type E2.2, 2.3 & E4.4). I30–32 (NM2) were associated with plant species occurring in meso-xeric calcareous meadows (Biscutella laevigata, Rhinanthus alectorolophus, Ornithogalum sp., Anthyllis vulneraria, Polygala vulgaris – EUNIS habitat type E1.26). At the upper margin of this group, I27 (NM1) could be interpreted as a wet variant of these subalpine meadows, with some species occurring in tall herb fringe habitats with waterlogging, i.e. forest/meadow ecotone in waterlogged conditions (Valeriana officinalis, Filipendula ulmaria, Bistorta officinalis - EUNIS habitat type E3.4) as well as with some underwood (Anemone nemorosa) and mesic meadow species (Lolium perenne, Luzula campestris). The last two NM1 sites (I21, I24) and one M site (I36) were associated with species typical of mesic meadows (e.g. Salvia pratensis, Arrhenatherum elatius, Leucanthemum vulgare, Trifolium pratense, Lotus corniculatus), but also with species from the herbaceous layer of mesic forests (Myosotis sylvatica, Primula elatior).

Fig. 3
figure 3

NMDS analysis, 2-D representation of the similarity (Jaccard’s index) among sites based on the occurrence of herbaceous plant species on sites where A. halleri was found. Black: metalliferous sites; red: NM1 sites; green: NM2 sites. Species scientific names (in grey) are abbreviated as in Table S2

Full size image

The last NM sites, I25 (NM1), I28, I29 and I33 (all NM2) were more dissimilar (more distant on the graphical representation) from each other and from the other NM sites, but still in the same part of the NMDS representation. These sites were associated with species found in ecotone habitats (Tephroseris longifolia: Janišová et al. 2012) (Fragaria vesca - EUNIS habitat E5: Woodland fringes and clearings and tall forb stands) in their mountain form (Senecio hercynicus, Aruncus dioicus, Valeriana tripteris), as well as rocky underwoods (Pseudoturritis turrita).

Polluted sites with A. halleri (M sites) were spread over the NMDS two-dimensional space (Fig. 3, in black). However, I36 was positioned within the “ecotone” group of NM sites, as described above. I12, I13, I15, I18 and I19 were relatively scattered over a large part of the NMDS 2-D space. Those M sites were mainly associated with species typical of mesic to eutrophic calcareous forest (e.g. meso to eutrophic forests G1.A Geranium robertianum, Aegopodium podagraria, Vinca minor, calcareous forests G1.72 and G1.63 Anemone hepatica, Cardamine heptaphylla, Mercurialis perennis, Polygonatum multiflorum, Ruscus aculeatus, Orchis mascula, Phyteuma spicatum) and of woodland fringes, clearings and tall forb stands (Solidago virgaurea, Hieracium lachenalii - EUNIS habitat E5) in a calcareous form (Dioscorea communis). I19 was positioned at the margin of this group, close to the last group of M sites, highlighting the mosaic of habitats of this site (Table 1), i.e. calcareous rocky grasslands and woods.

Finally, the last M sites, I16, I17 and I35 were highly distant from each other and from the other sites. They were all associated with xeric basophilous or calcareous grassland species (Globularia bisnagarica, Helianthemum nummularium, Stipa pennata, Sesleria caerulea, Teucrium chamaedrys – EUNIS habitat E1.27: sub-atlantic very dry calcareous grasslands), basophilous heaths (Erica carnea) and fringes (Cephalanthera longifolia, Vincetoxicum hirundinaria), as well as central European heavy-metal grasslands (Minuartia verna - EUNIS habitat E1.B3).

To summarize, A. halleri under non-metalliferous conditions (NM sites) therefore mostly occurred in mesic calcareous meadow sites, and in ecotones at the limit of forest and meadows. A. halleri in metal-polluted environments (M sites) occurred in more diverse and specific habitats, thriving over the whole vegetation succession going from calcareous xeric grasslands and calcareous mesic meadows, to calcareous forests and their fringes. NM1 and NM2 sites could barely be differentiated by this analysis, suggesting similar ecological habitats for A. halleri under non-polluted conditions in the two valleys.

Comparisons among groups for hyperaccumulation variables and abiotic parameters

In general, Zn and Cd shoot concentrations from our study (Table S1) were mostly among the lowest values obtained in A. halleri at the European scale (see Figure 1 from Stein et al. 2017). It was particularly true for shoot Cd concentrations, which were similar to those obtained by Stein et al. for populations from the Italian area (see Figure 4 from Stein et al. 2017). These authors suggested that these local values may characterize a particular genetic unit in this area, as shown by several authors (Pauwels et al. 2012; Šrámková et al. 2017). In our dataset, Zn and Cd shoot concentrations were highest in the M group, though shoot Cd concentrations did not significantly differ between M and NM1 groups (Fig. 4a and d). The latter result may rely on a few particular individuals in some NM1 populations (mostly individuals from I24 and I25, Table S1). Bioconcentration factors were significantly lower in M than in NM samples (Fig. 4c and f) since bioavailable soil Zn and Cd concentrations were higher in M than NM samples (Fig. 4b and e). There were no remarkable abnormalities of TME concentrations in non-metalliferous soils (except slightly elevated Cd concentrations in I27, see Table S1).

Fig. 4
figure 4

Boxplots and coefficients of variation for plant and soil variables related to Zn and Cd concentrations. Shoot Zn and Cd concentrations are provided in figures a) and d), respectively. Soil Zn and Cd concentrations are provided in figures b) and e), respectively. Bioconcentration factors (shoot metal concentration divided by soil metal concentration) of Zn and Cd are provided in figures c) and f), respectively. Coefficients of variation (%) are indicated within tables below each graph. M = metalliferous sites from “Val del Riso” valley; NM1 = non-metalliferous sites from « Val del Riso » valley; NM2 = non-metalliferous sites from « Sommaprada » or « Val Paisco » valleys. Letters above boxplots indicate significant differences at 5% significance threshold

Full size image

Regarding concentrations of other elements in soils, M group often significantly differed from NM groups, considering either all NM groups (Cu, Fe, K, Mg, Na, Pb: Fig. 5b, c, d, e, f and h) or only the NM2 group (Ca: Fig. 5a). Metalliferous soils were particularly impoverished in Fe and Na, but enriched in Ca, K, Cu and Pb (in addition to Zn and Cd: Fig. 4b and e). Element concentrations could also differ between NM1 and NM2 groups (Ca, Fe, Mg, P: Fig. 5a, c, e and g). Soils of NM2 sites were specifically enriched in Fe but showed low Ca and Mg concentrations.

Fig. 5
figure 5

Boxplots and coefficients of variation for soil variables related to other elements than Zn and Cd. Soil concentrations of macro-elements are provided in figures a), d), e), f) and g). Soil concentrations in micro-elements are provided in figures b), c) and h). Coefficients of variation (%) are indicated within tables below each graph. M = metalliferous sites from “Val del Riso” valley; NM1 = non-metalliferous sites from « Val del Riso » valley; NM2 = non-metalliferous sites from « Sommaprada » or « Val Paisco » valleys. Letters above boxplots indicate significant differences at 5% significance threshold

Full size image

Overall, the other environmental parameters showed less variability within (see coefficients of variation) and among edaphic groups than soil elemental concentrations (Fig. 6). There were no significant differences among groups for soil particle size distribution (Kruskal-Wallis test: 0.28 < p-values < 0.31), altitude (Kruskal-Wallis test: p-value = 0.057) and maximum solar irradiance (Kruskal-Wallis test: p-value = 0.23) (Fig. 6b, c, d, e, g). However, it is worth noting that there was a clear tendency for stonier soils and lower altitudes in Zn/Cd-rich habitats, and higher solar irradiance in NM2. Only pH and absolute humidity values were significantly contrasted among edaphic groups, with mostly lower values in NM2 (Fig. 6a and f). The decrease of pH values for NM2 was easily explained by the occurrence of a pine forest among surveyed populations (I33, see Table 1).

Fig. 6
figure 6

Boxplots and coefficients of variation for other soil variables (a, b, c and d), topographic (e) and climatic variables (f, g). Coefficients of variation (%) are indicated within tables below each graph. M = metalliferous sites from “Val del Riso” valley; NM1 = non-metalliferous sites from « Val del Riso » valley; NM2 = non-metalliferous sites from « Sommaprada » or « Val Paisco » valleys. Letters above boxplots indicate significant differences at 5% significance threshold

Full size image

PLSR on hyperaccumulation variables

The four regression analyses performed on the whole data set revealed that soil Zn and Cd concentrations as well as soil particle size distribution always significantly affected the four hyperaccumulation variables (Table 2). Other soil elements (like Ca, Na or Pb), soil pH and absolute humidity were also involved for some of the four hyperaccumulation variables.

Table 2 Results of PLS regression on standardized data
Full size table

The 12 regression analyses performed on the different groups showed contrasted results among the three groups and among the four hyperaccumulation variables (Table 2). Considering variables related to Cd hyperaccumulation ([Cd] and BF_Cd), the most striking result was the incapacity of the statistical procedure to extract explaining factors in NM1 and NM2 groups. For the M group, soil Cd concentration showed positive effect on Cd hyperaccumulation values but negative effect on the bioconcentration factor. Several other elements showed a significant effect on Cd hyperaccumulation (Cu, K, Pb, Zn) whereas only soil Pb concentrations acted on the bioconcentration factor. In addition, soil particle size distribution seemed to explain Cd hyperaccumulation variation. Topographic and climatic variables also influenced the bioconcentration factor of Cd (Table 2).

Considering Zn hyperaccumulation variables ([Zn] and BF_ Zn), all the abiotic parameters significantly affected shoot concentrations or bioconcentration factors in at least one group (Table 2). However, since they were not based on the same numbers of components, the different regression models were not fully comparable. For this reason, comparisons among models only focused on the significance of parameters in each model (based on VIP values) and signs of regression coefficients. The effects of the different parameters could therefore be specified as “0” for no effect, “+” for positive significant effect, “-” for negative significant effect (Table 3). Consequently, the regression models for [Zn] and BF_ Zn could be compared and opposed the three groups of populations in different ways, depending on the distribution of effects (0, +, −) among abiotic parameters: M = NM1 / NM2, NM1 = NM2 / M, M = NM2 / NM1, all groups different, all groups similar (Table 3). Significant opposite effects in different groups (e.g. +/− or −/+) were rarely observed, whereas positive or negative significant effect in one or two groups and no significant effect in the other(s) (e.g. 0/+ or 0/−) was a frequent situation. Interestingly, for shoot Zn concentration ([Zn]), soil Cd concentration clearly opposed the two valleys (+/− for M = NM1 vs NM2), whereas soil Zn concentration affected all groups of populations in a positive way (+ for M = NM1 = NM2) (Table 3A). For bioconcentration factor of Zn (BF_Zn), soil Zn concentration also opposed the two valleys (−/+ for M = NM1 vs NM2), while all the other parameters did not have so contrasted effects (always one “0” effect in the group comparison) (Table 3B).

Table 3 Summary of PLS regression results showing the relative effect of significant abiotic parameters on groups of populations
Full size table

Discussion

Habitat heterogeneity among and within M and NM groups

In this study, the habitat of A. halleri, a pseudometallophyte model plant, was described in several sampling sites varying for many abiotic environmental parameters. The studied region was of particular interest for several reasons. First, this region includes former mining and industrial sites in a landscape where A. halleri naturally occurs. This made it more likely to observe metalliferous and non-metalliferous sites in close geographical proximity, so that corresponding populations may share common recent history (Pauwels et al. 2012). In contrast, previous studies investigating the ecogeographic distribution of metal-related traits involved metalliferous sites in the margin of the species range (if not outside), at a certain geographic distance of any non-metalliferous sites (Meyer et al. 2010; Pauwels et al. 2006). In such situation, genetic relationships between metallicolous and non-metallicolous populations may also result from complex demographic histories, so that phenotypic divergence among populations may involve non-adaptive processes (Pauwels et al. 2012). Second, some populations from the region we investigated have been recently mentioned in the literature (Meyer et al. 2015; Stein et al. 2017), but the local ecological context has never been reported anywhere.

In 23 populations of this region, we investigated not only soil elemental concentrations, but also several abiotic parameters (soil pH and structure, climate, topography), in order to detect spatial discrepancies at a local scale.

In order to better interpret results, we gathered sampling sites into three groups, depending on their geographic location (M-NM1 vs NM2) and their close proximity to a mining or industrial area (M vs NM1-NM2). Analyses of soil metal concentrations further supported the occurrence of true metalliferous soils on metallurgic sites (M), whereas other sites were largely non-metalliferous (NM1 and NM2) (see Table S1).

The ecological description of metalliferous and non-metalliferous habitats provided further insights into our understanding of the ecology of A. halleri. As far as we know, metalliferous environments were generally described as exhibiting harsh living conditions, where TME concentrations are very toxic, and where soils are weakly structured, free draining and moisture deficient, and show low nutrient status (Macnair 1987, 1997). On the contrary, non-metalliferous habitats of pseudometallophytes have been little described, except for N. caerulescens populations from Grand Duchy of Luxemburg (Molitor et al. 2005). As expected, we observed high TME concentrations in all sampling sites close to metallurgic activities. We also observed lower Fe and Na concentrations in metalliferous soils, without being able to confirm a deficiency of these elements (Fig. 5c, f). Metalliferous soils also tended to be stonier than non-polluted soils from the same region (Fig. 6d). However, contrary to what is currently admitted, in our study metalliferous soils were not systematically impoverished in essential elements, since Ca, K and P concentrations were higher or not so different from those in non-metalliferous soils (Fig. 5a, d, g). We also evidenced, from both floristic and abiotic data (Figs. 2 and 3), a stronger heterogeneity among metalliferous habitats of A. halleri than among non-metalliferous habitats. Under polluted environmental conditions, A. halleri actually thrives in a variety of habitats covering all stages of the succession, i.e. from calcareous grasslands to forests. In addition, M sites were characterized by a large variability of altitudes and solar irradiance values (Fig. 6e and g). In comparison, non-metalliferous habitats showed higher homogeneity. Under non-polluted environmental conditions, A. halleri was indeed mostly found in calcareous mesic meadows, as well as in ecotone habitats at the fringe of forests and meadows. Notably, non-metalliferous habitats were only slightly differentiated depending on their valley (i.e. NM1 vs NM2), since only 15.17% of variance was explained by the NM1-NM2 separation in PCA analysis (Fig. 2a). This differentiation was due to some abiotic parameters (Mg and Fe soil concentrations, absolute humidity, solar irradiance), but was not confirmed by the study of plant communities (Fig. 3), showing that abiotic variations have a rather low impact. The fact of being more heterogeneous is therefore a peculiarity of metalliferous habitats. This has already been suggested for Noccaea caerulescens populations but only on the basis of mineral soil concentrations and by comparing four populations (Dechamps et al. 2011). Moreover, such heterogeneity within and among metalliferous habitats suggests that the current and binary classification into metalliferous and non-metalliferous, which is only based on TME soil concentrations, may be misleading.

Ecological niche extension during M sites colonization

The range of values for a set of abiotic and biotic environmental parameters which can be reported from a species geographical range is assumed to reveal its realized ecological niche (Pironon et al. 2017). Accordingly, the colonisation of new habitats showing environmental conditions outside the original species range can be related to evolution of adaptive ecological traits and/or ecological niche extension (Holt 2003). In A. halleri, the ecological differences among metalliferous and non-metalliferous sites, as well as the higher heterogeneity among metalliferous sites, suggest that their colonisation by A. halleri might have significantly extended the realized ecological niche of the species, in several environmental dimensions. The extension of the ecological niche during the colonisation of industrial sites was already suggested for A. halleri in Northern France, where the species only occurs in polluted areas and was locally described as non-native (Berton 1946; Van Rossum et al. 2004). In particular, ecological niche modelling from topographic (1) and climatic (18) data previously suggested that populations from Northern France settled outside of the natural climatic niche of the species (Fig. S4, Pauwels et al. 2012), suggesting population maintenance despite unsuited climatic conditions. The reason why metal pollution could facilitate the naturalization of the species outside its natural niche is still elusive. Through interactions with other abiotic or biotic environmental parameters, TME may reduce stress in otherwise unsuitable conditions and allow tolerant plants to maintain sufficient performance for population maintenance. For example, TME pollution may modify plant communities and reduce interspecific competition and/or herbivory. Thus, it was demonstrated that herbivores’ densities were reduced in metalliferous environments, to such extent that metallicolous populations evolved towards constitutively lower glucosinolate concentrations (Noret et al. 2007). Alternatively, adaptation to TME pollution may involve concomitant adaptation to other abiotic environmental constraints, like, for example, drought that can be caused by the stony structure of many metalliferous soils. Such tolerance to multiple abiotic stresses was demonstrated for some halophyte species (for review, see: Ben Hamed et al. 2013). Subsequently, adaptation to ecological conditions that were initially outside the species niche may allow species range extension (Holt 2003). Testing such a hypothesis will require assessing the tolerance of plants to various stress conditions in controlled conditions. A mean to investigate more precisely to which ecological (biotic or abiotic) constraints plants may respond during the colonization of new habitats would be the study of the variation in leaf functional traits along a metal contamination gradient, as it was performed on serpentine and cupriferous soils (Kazakou et al. 2008; Lange et al. 2017). It has been indeed suggested that functional trait-based approaches may theoretically inform about the ecological drivers affecting plant functional responses (Enquist et al. 2015). These approaches can also be extended to the whole plant community, by measuring both leaf functional traits and species abundances, as it was achieved along copper-cobalt gradients (Delhaye et al. 2016). What is known so far for A. halleri, is that this species was somehow preadapted to metal-pollution through ancient and constitutive overexpression of a metal transporter conferring high tolerance to Zn and Cd (the Heavy Metal ATPase 4: Meyer et al. 2016). Hence, the species could have taken advantage of this preadaptation to colonize recently created and, for any reason, favourable habitats.

Potential selective pressures acting on Zn and Cd hyperaccumulation

Phenotypic and environmental data confirmed that hyperaccumulation abilities were mostly structured among M and NM sites, although we could not statistically distinguish NM1 sites from M ones for shoot Cd concentrations. Our results are congruent with previous studies showing contrasting hyperaccumulation capacities of metallicolous and non-metallicolous individuals in situ and/or in controlled conditions under Zn treatment (Bert et al. 2002; Stein et al. 2017). In these studies, Zn hyperaccumulation capacities of metallicolous individuals were mostly inferior to that of non-metallicolous individuals under Zn treatment.

Then, we assumed that PLSR analyses could reveal whether a particular combination of abiotic parameters could significantly affect these Zn or Cd hyperaccumulation abilities. PLSR results were generally contrasted between Zn/Cd shoot concentration and bioconcentration factors (Table 2). Since bioconcentration factors should represent hyperaccumulation abilities independently from soil concentration values, we considered that particular attention has to be paid to their response patterns in PLSR analyses.

Overall, our results on the whole dataset suggested that the abiotic parameters that explained phenotypic differences for all hyperaccumulation variables were also those that mostly distinguish metalliferous and non-metalliferous habitats, i.e. soil Zn and Cd metal concentrations, and soil particle size distribution as a correlated parameter (Table 2 and Fig. 2c). Considering group-specific results, we first observed that Cd hyperaccumulation variables were badly predicted by PLSR models in non-metalliferous habitats. In a comprehensive analysis of 165 A. halleri populations, Stein et al. (2017) also showed that bioavailable soil Cd concentrations explained only a small part of Cd hyperaccumulation, in particular when soil Cd concentrations were low. Such a result suggests several hypotheses. It may imply that the evolution of Cd hyperaccumulation could be neutral in non-polluted habitats. Alternatively, it may imply that, in the absence of pollution, the evolution of Cd hyperaccumulation is not strongly affected by the abiotic parameters we recorded. Several authors actually proposed that Cd hyperaccumulation could have evolved under biotic selective pressures, more precisely as a defence against herbivores (Kazemi-Dinan et al. 2014; Plaza et al. 2015). However, this evolutionary hypothesis should be more plausible in non-metalliferous habitats than in metalliferous ones where herbivores only occur in very low densities (Dechamps et al. 2008; Noret et al. 2007). In metalliferous habitats, our results conversely suggest that [Cd] in shoots and BF_Cd can be predicted by soil [Cd] and soil [Pb] (alone for BF_Cd, together with other parameters for shoot [Cd]), two pollutants that generally co-occur in calamine sites. This would suggest that when soil [Cd] is high, i.e. when Cd toxicity is elevated, it may represent a more powerful driver of the evolution of Cd-related traits that when it is low.

In our study, variation in shoot Zn concentrations or bioconcentration factors of Zn were quite well predicted by PLSR models in all groups of populations. This suggests that the abiotic parameters recorded here may partly structure the polymorphism of Zn hyperaccumulation in natural populations. Interestingly, this is true either in metalliferous or non-metalliferous habitats, suggesting that among-habitat variation in abiotic parameters may influence shoot [Zn] and BF_Zn whatever the level of soil metal pollution. Soil Zn concentration was the only parameter that showed VIP >0.8 in the six regression models. This confirmed that it may play a special role in structuring Zn hyperaccumulation capacities. This parameter positively affected shoot Zn concentration in all groups of sites (Table 3A), probably because Zn uptake mainly tends to naturally increase with Zn availability in soil (Stein et al. 2017). However, soil Zn concentration negatively affected BF_Zn in M-NM1 sites, while M and NM1 plants displayed significantly different Zn hyperaccumulation capacities (Fig. 4a, c), but positively affected BF_Zn in NM2 sites (Table 3B). Such a discrepancy between M-NM1 vs NM2 individuals was actually observed for several abiotic parameters (Table 3A, B). The reason of this phenomenon is unclear. At least, it would be tempting to assume that geographic distances should impact gene flow, and that higher gene flow between M and NM1 sites within Val del Riso valley could influence among-population genetic and phenotypic differentiation. As an example, significant influence of neutral genetic structure on PLSR analyses, performed to explain phenological data by ecological variables, was evidenced in Arabidopsis thaliana (Brachi et al. 2013).

It is also worth to indicate that, for bioconcentration factors of Zn, results opposing M and NM groups were significant for two abiotic parameters. First, soil Cd concentration had a negative influence on bioconcentration factors on non-polluted sites (NM1 and NM2) comparing to polluted sites (Table 3B). Competition between Zn and Cd and the subsequent decrease of Zn accumulation in the presence of Cd in the growth medium was already observed in controlled conditions by several authors (Assunção et al. 2003; Escarré et al. 2013; Roosens et al. 2003). Absolute humidity, which was not correlated to soil Cd concentration (Fig. 2c), also seemed to have a positive effect on bioconcentration factor of Zn only in non-metalliferous habitats (Table 2). For some abiotic parameters, we could have expected a significant effect according to the description of their role in literature. In particular, several studies highlighted the influence of pH on either Zn and Cd availability or Zn and Cd absorption (see for example: Degryse et al. 2009; Maxted et al. 2007). In our study, pH had a significant negative effect on shoot Zn concentration in plants from the M group, which is in accordance with the decreasing of metal availability when pH increases (Degryse et al. 2009). Anyway, it is important to note that PLSR analyses do not provide evidence of a causal effect of explaining parameters on a dependent variable, so that interpretation of results should always be cautious.

An interesting perspective to these PLSR analyses could be to perform Generalized Additive Model (GAM) analyses on the more informative abiotic parameters evidenced here. Such approach should explain how Zn and Cd accumulation values, or Zn and Cd bioconcentration factors, may change with A. halleri niche parameters, and may potentially reveal threshold effects.

Conclusion

In the present study, using both environmental and floristic data, we clearly demonstrated that metalliferous habitats are not just metal-polluted but also more ecologically diverse than non-metalliferous ones. The colonization of metalliferous habitats may thus correspond to an enlargement of the ecological niche of A. halleri, potentially allowed by its tolerance and hyperaccumulation species-wide properties. We could not exclude that soil Zn and Cd concentrations may represent selective agents, acting on Zn hyperaccumulation in particular. However, other abiotic parameters may also shape the natural variation of this trait. Combinations of abiotic parameters, specific to each habitat (M, NM1 or NM2), probably represent distinct regimes of selective pressures. Our results pave the way to new assumptions about potential selective pressures for the evolution of metal accumulation. This approach now needs to be completed by direct experiments in controlled conditions. For instance, experimental evolution would be a convincing way to prove the action of any abiotic parameters as selective pressure for the evolution of metal-related traits.