Patterns of genetic divergence among populations of the pseudometallophyte Biscutella laevigata from southern Poland

Pseudometallophytes are model organisms for adaptation and population differentiation because they persist in contrasting edaphic conditions of metalliferous and non-metalliferous habitats. We examine patterns of genetic divergence and local adaptation of Biscutella laevigata to assess historical and evolutionary processes shaping its genetic structure. We sampled all known populations of B. laevigata in Poland and analyzed respective soil metal concentrations. For genotyping we used nine nuclear microsatellite loci. Population genetic pools were identified (Bayesian clustering) and we estimated genetic parameters and demographic divergence between metallicolous and non-metallicolous populations (ABC-approach). Populations clustered into two groups which corresponded to their edaphic origin and diverged 1,200 generations ago. We detected a significant decrease in genetic diversity and evidence for a recent bottleneck in metallicolous populations. Genetic structure was unrelated to site distribution but is rather influenced by environmental conditions (i.e. soil metal concentration). The intriguing disjunctive distribution of B. laevigata in Poland results from a fragmentation of the species range during the Holocene, rather than recent long-distance-dispersal events. The genetic structure of populations, however, continues to be modified by microevolutionary processes at anthropogenic sites. These clear divergence patterns promote B. laevigata as a model species for plant adaptation to polluted environments.


Introduction
The adaptation of organisms to natural and/or anthropogenic changes in their physical environment has raised considerable interest among plant ecologists and evolutionary biologists. Special attention has been given to living organisms occurring on sites with high concentrations of Metal Trace Elements (hereafter called "metals"). Most organisms are unable to cope with such harsh environmental conditions, but some plant species manage to colonize, develop, and reproduce in metal polluted habitats. Those species are thus defined as metal tolerant (Macnair 1987). Metal tolerance may have evolved on either naturally metalenriched soils or on anthropogenic sites that have been generated by waste deposits from ore mining or metal smelter activity (Ernst et al. 1990;Peralta-Videa et al. 2009). In contrast to natural metalliferous soils, industrial sites are younger and generally show higher metal concentrations. These differences in site history may affect the evolutionary processes and selection pressure dynamics. Adaptation is expected to evolve faster on artificially polluted sites, but with large differences observed between taxa (Ernst 2006). Hence, the evolution of metal tolerance on anthropogenic contaminated soils has been extensively studied for more than half a century to better understand plant adaptation to restrictive environments over ecological times (Antonovics et al. 1971;Linhart and Grant 1996).
Plants occurring on metal contaminated soils usually belong to pseudometallophyte species developing both metallicolous (M, on metal-enriched soils) and nonmetallicolous populations (NM, on non-polluted soils; Lambinon and Auquier 1964). Such contrasting edaphic types are expected to be genetically distinct (Macnair et al. 2000). Firstly, local adaptation to metal-polluted soils is expected to result in a higher frequency of tolerant genotypes in M compared to NM populations. Accordingly, several phenotyping experiments reported higher mean metal tolerance levels in M compared to NM populations (e.g. Assunção et al. 2003;Meyer et al. 2010;Pauwels et al. 2006). Secondly, considering that M populations of pseudometallophyte species occurring on anthropogenic polluted sites result from recent colonization and adaptation events, their level of neutral genetic diversity is expected to be reduced due to selection, founder and bottleneck effects during the colonization of metalliferous soils (Lefèbvre and Vernet 1990;Wu 1990). This rationale has been confirmed in several studies (Deng et al. 2007;Mengoni et al. 2001;Meyer et al. 2010;Nordal et al. 1999;Wojcik et al. 2013). Other investigations, however, found comparable levels of neutral genetic variability in both M and NM populations (Baumbach and Hellwig 2007;Bizoux et al. 2008;Pauwels et al. 2005), or even slightly higher level within the M group (Słomka et al. 2011). Such contrasting patterns suggest that a simple comparison of edaphic types may not be sufficient to understand the structure of genetic diversity within and among populations of pseudometallophyte species. Additional factors that could influence the overall level of genetic diversity should be considered, including i) the age of polluted sites, ii) the level of soil metal concentration and associated selective pressure, iii) the level of genetic connectivity between M and NM populations, iv) the possibility of successive colonization events from adjacent NM populations, and v) some species-specific features such as clonal growth, sensitivity to abiotic stress, and constitutive (i.e. species wide) or population-specific metal tolerance (Bickham et al. 2000;Kang and Chung 2000;Vekemans and Hardy 2004). So far it has remained largely unclear whether historical, ecological or biogeographic factors could influence the pattern of genetic structure among M and NM populations of pseudometallophytes (Ye et al. 2012).
Southern Poland with its mining sites and abundance of metal tolerant species is an ideal location to study factors influencing pseudometallophyte population genetics (Zając 1996). Mining and processing of rich zinclead (Zn-Pb) ore in the Bolesław -Olkusz region (Cracow -Silesian Upland) dating back to the 13th century has created some of the largest and most polluted anthropogenic metalliferous sites in Europe, with waste heaps and dust deposits of different age, composition and metal concentrations (Cabala et al. 2009). A specific calamine flora (Grodzinska et al. 2000) has evolved in response to the selection pressure from elevated metal concentrations (Zn, Pb, Cd) in soils. The majority of calamine species are native and widespread in the area surrounding the heaps (Jędrzejczyk-Korycińska 2006). For those species, M populations can reasonably be assumed to derive from neighboring NM populations (Grodzinska et al. 2000). However, calamine flora also includes mountain species that in lowland areas occur almost exclusively on polluted soils (Grodzinska et al. 2000;Zając 1996). As a result, the distribution of such species is highly disjunctive and the closest NM populations are located~100 km away in the Tatra Mts. Under such circumstances, lowland M populations may either result from long distance dispersal (LDD) events postdating pollution or represent relic populations from a formerly larger distribution range prior to the Last Glacial Maximum (LGM). Either way, the genetic exchange between M and NM populations must be strongly limited.
Biscutella laevigata L. (Brassicaceae) is a typical but poorly investigated pseudometallophyte, which is present in the metal-polluted habitats of southern Poland (Wóycicki 1913). Its ability to tolerate large amounts of some metals in foliar tissues (Pošćić et al. 2013) shows an interesting potential as a model species to study the adaptive colonization of metalliferous sites. This adaptive potential and tolerance to stressing ecological conditions were previously demonstrated for serpentine alpine populations (Gasser 1986). Biscutella laevigata is a perennial, strictly outcrossing species (Olowokudejo and Heywood 1984) which is widespread across Central and Western Europe and reaches its northern distribution limit in Poland (Jalas et al. 1996;Tremetsberger et al. 2002). There, B. laevigata has a clearly disjunctive distribution. It is mostly present on non-metalliferous sites in the Western Tatra Mts. However, it also grows on calamine waste heaps (i.e. metalliferous sites) in the Bolesław -Olkusz region, located~100 km north of this mountain range. Furthermore, one occurrence on a nonmetalliferous lowland site in Zagórzyce (Przemyski and Piwowarczyk 2012; subsequently referred to as PL10) has recently been mentioned. Populations from both the lowland waste heap sites and the Tatra Mts have been identified as diploids (2n=18) in a caryological study (Skalińska 1950) and varied in morphological and physiological traits, as well as metal tolerance levels (Wierzbicka and Pielichowska 2004). It has remained unclear, however, whether lowland populations result from relatively recent LDD events or if they are glacial relic.
In the present study, we investigated the genetic variability of B. laevigata among all 16 known low and high elevation provenances from locations in southern Poland, where highly variable metal concentrations in soil were reported (Szarek-Łukaszewska and Niklińska 2002;Babst-Kostecka et al. in review). We surveyed the nuclear genome diversity to examine population-level microevolutionary processes in the context of environmental and geographic variation. The specific aims of this study were to i) assess the genetic variability of B. laevigata, ii) compare patterns of genetic diversity among populations growing on recent (anthropogenic) and ancient (natural) habitats to test for a reduction in genetic variability after colonization of polluted areas, and iii) evaluate a demographic scenario and the timing of divergence between M and NM populations. By addressing these points, our study provides new insight into the origin and biogeographic characteristics of M and NM populations of pseudometallophyte species.

Sampling and chemical analysis
Leaf material from 522 diploid individuals was sampled every 3 m along 100 m long transects at all of the 16 known (sub-)alpine locations of B. laevigata in the Tatra Mts (i.e. sites PL8, PL9, PL11-PL13, SK14, SK15, PL16) and their northern foothills (PL1-PL7, PL10; Fig. 1 and Table 1). Locations were categorized as Fig. 1 Location of study sites. Open and closed squares represent sites from non-industrial and industrial areas metalliferous and non-metalliferous according to the total Zn concentration in soil. Thereby, we considered 300 mg . kg −1 as the threshold concentration (Bert et al. 2002). Sample size varied from 28 to 37 individuals per location and distances between sites ranged from 1 km to more than 150 km. Rhizospheric soil (5 samples per site) was sampled to a depth of 10 cm using a cylinder of 7 cm diameter. Soil pH was measured in distilled water by means of the potentiometric method (Wąchalewski 1999). Total concentrations of Zn, Pb and cadmium (Cd) were determined using atomic absorption spectrometry. The extractable fraction was determined by vortexing soil solution with EDTA (Escarré et al. 2000). Differences between average total concentrations of Zn, Pb and Cd, and their extractable fractions at metalliferous vs non-metalliferous sites were tested by means of Student's t-test. The Bartlett test was performed to compare the homogeneity of variances between natural and anthropogenic sites.

Molecular analysis
The collected leaf material was dried for 24 h at 55°C prior to DNA extraction from 10 to 15 mg of dry material using the NucleoSpin® 8/96 Plant kit (Macherey-Nagel®). PCR amplification was performed using 1:100 dilutions.

Analyses of genetic structure and diversity
The mean allelic richness (A S ) within each site was estimated over loci using FSTAT 2.9.3 (Goudet 2001). The rarefaction method (Kalinowski 2004) was used to standardize estimates to the smallest sample size (n=28) collected from population PL1. The expected (H E ) and observed (H O ) heterozygosity as well as Wright's inbreeding coefficient (F IS ) were estimated for each site using Genetix 4.05 (Laboratoire Génome et Populations, CNRS UPR 9060, Université de Montpellier II, Montpellier, France). Differences  (Goudet 2005). Differences between M and NM groups were tested using randomization procedures (1,000 permutations of populations between edaphic types). We estimated the number of homogeneous gene pools (K) that best explained the pattern of genetic structure at microsatellite loci using Bayesian modelbased clustering implemented in Structure 2.3.2 (Pritchard et al. 2000). Analyses were carried out under the admixture model, assuming an identical parameter of individual admixture across clusters and a uniform prior. For each K value within the range of 1-16 we performed 10 independent runs with 100,000 iterations after a burn-in period of 10,000 generations, with no prior information on the origin of individuals. Finally, from the mean likelihood log e (K) of replicates we computed the ΔK statistics to identify the most likely number of clusters (Evanno et al. 2005). Additionally, a Neighbor-Joining (NJ) population tree was drawn from Cavalli-Sforza and Edwards chord distances, based on allele frequencies and using the TreeView software (http://taxonomy.zoology.gla.ac.uk/rod/treeview.html). The node robustness was evaluated using 10,000 bootstrapped replicates.
The contribution of geographic distances to differences in the genetic structure (isolation by distance patterns, IBD; Wright 1949) between populations was assessed by correlating the matrices of log-transformed geographic distances between pairs of populations and linearized pairwise genetic distances (F ST /(1-F ST )) using a Mantel test (Mantel 1967). Pairwise F ST values were calculated using 1,000 random permutations of individuals in Genetix 4.05 (Laboratoire Génome et Populations, CNRS UPR 9060, Université de Montpellier II, Montpellier, France). Separate Mantel tests with 100,000 random permutations each were performed for the entire data set, as well as for the M and NM populations separately, using the R-package ade4 (Dray and Dufour 2007).

Approximate Bayesian Computation (ABC) analysis
A demographic scenario of divergence between M and NM populations was modeled and the simulated datasets were compared to observed data by means of an ABC approach (Beaumont et al. 2002). Three parameters were considered to describe the demographic model: effective population size of M (N M ), effective population size of NM (N NM ), and divergence time in generations (t G ). Using the program DiyABC v2.0 (Cornuet et al. 2008), we simulated 500,000 multilocus datasets based on a demographic history that describes the model. For each simulation, the parameter values were drawn from their prior distributions. The uniform prior distributions were defined over the intervals 10 -200,000 for N M and N NM , and 10 -5,000 for t G . The mutation rate and the model for microsatellites were set to default. We generated distributions of 11 summary statistics of genetic diversity (mean number of alleles per locus: A M , A M ; mean expected heterozygosity: H eM , H eNM ; mean allele size variance: V M , V NM ; F ST between pairs of population samples, classification index (also called mean individual assignment likelihood): LIK M , LIK NM ; shared allele distance: DAS; (dμ) 2 distance between two samples) by coalescent-based simulations. The posterior distributions of the parameters for the model were estimated on the 2,000 replicate simulations with the smallest associated Euclidean distance (δ=║ S obs -S sim ║) using a non-linear regression after the application of a logit transformation on parameter values (Beaumont et al. 2002). Finally, following Cornuet et al. (2010), we employed the model-checking function of DiyABC to assess the goodness-of-fit between each model parameter -posterior combination and the observed dataset. Computations based on 1,000 simulated datasets and 13 summary statistics (i.e. A, H e , V, mean M GW index (Garza-Williamson's M) across loci, LIK, DAS, (dμ) 2 distance), generated a posterior cumulative distribution function for each statistic together with the estimated probability (P values) of these distributions to fit the original observations. A principal component analysis (PCA) was performed in the space of summary statistics, considering 5,000 datasets simulated with parameter values drawn from the prior. The observed dataset, as well as the 1,000 datasets simulated from the posterior distributions of parameters were passively projected into the PCA plane.

Detection of recent bottlenecks
We tested for a recent occurrence of bottlenecks in each of the 16 populations using Bottleneck 1.2.2 ) under a two-phase mutation model (TPM; Di Rienzo et al. 1994) and using 1,000 iterations. Three different percentages of mutations that follow a strict stepwise mutational (SMM) process (P SMM =5 %, 30 % and 60 %) and five values of the variance in size of multistep mutations (V TPM =2, 4, 9, 16, and 25) were tested. The V TPM were specified on the basis of the observed ranges of microsatellite allele sizes detected in our sampling (data not shown) and corresponded to average multistate mutational jumps of~1, 2, 3, 4, or 5 steps (Di Rienzo et al. 1994). We used the following two criteria to detect bottlenecks: 1) significant deviation from mutation drift equilibrium (H 0 =none-bottleneck; p<0.05, sign test; Luikart and Cornuet 1998) and 2) high probability of heterozygosity excess (p >0.5, Wilcoxon sign-rank test; .

Microsatellite diversity and population structure
All nine nuclear microsatellite loci investigated in this study were highly polymorphic, showing between 4 and 25 alleles per locus. Mean allelic richness (A S ), observed (H O ) and expected heterozygosity (H E ) ranged from 4.506 to 8.872, 0.284 to 0.658 and 0.309 to 0.743, respectively, and were significantly lower in M vs NM populations (p<0.001, t-test; Table 2). The inbreeding coefficient (F IS ) values were low, as expected for a selfincompatible species, but significantly higher than zero in all populations, which could relate to local biparental inbreeding. The level of genetic differentiation varied between populations, with a general tendency toward lower overall differentiation in M than in NM populations (p<0.001, t-test).

Bayesian clustering, IBD patterns and ABC approach
Based on the highest modal value of the ΔK statistics, the Bayesian clustering analysis revealed two sets of populations (Fig. S1). Considering K=2 as the real number of clusters, we observed a clear separation between populations from the metalliferous (cluster 1, C1) and non-metalliferous sites (cluster 2, C2; Fig. 2a). At K=2, distribution values (Q Ji ) were generally high and individuals from a given location were consistently assigned to a common gene pool. In some NM populations (in particular PL10, PL13-PL16), however, Q Ji values were lower and heterogeneous for several individuals (samples were not clearly assigned to the cluster C2), suggesting admixture. The NJ tree further revealed the existence of two distinct genetic groups (Fig. 2b). The inferred tree topologies were robust with most bootstrapped values greater 70 %. We did not find significant patterns of isolation by distance (IBD) for pairwise comparisons among all investigated populations (Fig. S2, r=−0.126, p=0.961, Mantel test), suggesting no contribution of geographic distance to genetic structure among edaphic types. Similarly, no correlation between genetic and geographic distances was observed for pairwise comparisons among samples within the M group (r=−0.288, p= 0.950). Within the entire NM group, the correlation was positive and significant (r=0.881, p=0.021) but turned non-significant if the only lowland location PL10 was removed from the analysis (r=0.205, p= 0.236). Overall, this suggests that gene flow within clusters is not restricted by geographic distances among sites.
Assuming a simple scenario of demographic divergence between M and NM populations, we inferred the population history and dated the median time of divergence between M and NM populations (t G ) back to 1,200 generations ago (Table 3). The median value of the mutation rate of SSR and SNI at the examined loci was estimated at 2.93×10 −4 and 4.56×10 −8 , respectively. Differences between the prior and posterior density curves for all demographic parameters were pronounced Metallicolous populations are marked in bold, non-metallicolous in normal font. C1 and C2 are the clusters of studied populations, according to their genetic distances with clearly peaked posteriors (Fig. S3), indicating that the genetic data were informative for parameters under the proposed scenario. Observed values of the summary statistics did not significantly differ from the simulated values based on parameter values drawn from the posterior distributions for the tested scenario (data not shown). In the PCA, the observed data matched the values obtained from the posterior and prior distributions (data not shown), thus further indicating a good fit of the posterior distribution based on the tested scenario.

Bottleneck detection
All M populations except PL2 consistently indicated a recent bottleneck event (Table S2). By contrast, patterns detected in NM populations changed, depending on the implemented mutational models (varying P SMM and V TPM values). In models with the lowest V TPM values (corresponding to average multistate mutational jumps of~1, 2 and 3), some NM populations presented a pattern coherent with a mutation-drift equilibrium (PL7, PL8, PL13, PL16), while others exhibited a bottleneck signature (PL9, PL12, SK14 and SK15). If the highest V TPM values were considered, some NM populations (PL7, PL8, PL13, PL16) showed high probability of heterozygosity deficiency (P(Hd)>0.7, Wilcoxon sign-rank test) suggesting a recent population expansion or the introduction of new alleles by immigrants.

Discussion
The present study highlights patterns of genetic diversity and population structure of all known M and NM populations of B. laevigata in Poland using nuclear microsatellite markers. This investigation showed that populations clustered into two groups which mainly reflected their edaphic origin. This remains true even considering the only lowland NM population (PL10) that was assigned to the geographically distant NM cluster of alpine sites.
In case of metallophytes, influences of edaphic conditions and soil metal contamination have often been shown to exceed the impact of geographic distance on genetic structure (Deng et al. 2007). Based on our results, similar trends can be proposed for B. laevigata, although this may not be fully conclusive in the absence of other lowland NM populations. Interestingly, an earlier investigation on a continuous tetraploid population of B. laevigata located in the Alps (Parisod and Christin 2008) similarly suggested that environmental rather than geographical conditions played a major role in shaping the genetic structure between alpine and subalpine habitats.
In contrast to multiple previous studies on pseudometallophytes, which failed to detect evidence for strong founder events in M populations (Pauwels et al. 2005 and references therein), our investigation suggested that the presence of genotypes in metal-contaminated environments is associated with genetic features supporting a founder effect in B. laevigata. Firstly, the genetic diversity observed in M populations was significantly lower compared to NM populations. Secondly, most M populations showed a strong and consistent bottleneck signal that was retained at all loci, suggesting recent demographic events . Considering that anthropogenic habitats are also of recent origin, genetic results strongly suggest that the occurrence of B. laevigata in lowland locations in the Bolesław -Olkusz region is associated with the colonization of metal-polluted sites. Moreover, the increased metal tolerance levels observed in M populations (Wierzbicka and Pielichowska 2004) also suggest that those populations experienced local adaptation to metal exposure, a phenomenon that could have further decreased genetic diversity. In summary, it appears that colonization of metal-polluted areas has been a dominant evolutionary driver of B. laevigata populations in southern Poland and has significantly modified the genetic structure of the species in recent times. A recent colonization of the Bolesław -Olkusz region would suggest LDD event from the Tatra Mts. In addition, we did not observe any evidence for a general pattern of IBD (i.e. for the balance between gene flow and genetic drift; Hutchison and Templeton 1999), which is likely a consequence of the uneven habitat distribution in southern Poland. However, B. laevigata has no biological capacity of LDD of seeds (Parisod and Bonvin 2008), in particular considering that the investigated edaphic types are located more than 100 km apart. Still, we cannot fully exclude the occurrence of LDD in the species history, as some molecular studies have shown patterns of highly disjunctive distributions via LDD, even for taxa characterized by limited seed dispersal capacity (e.g. Escudero et al. 2010;Schönswetter et al. 2004;Tribsch et al. 2002). Moreover, LDD could have been favored by human activities. Indeed, as B. laevigata grows on metal-enriched sites in lowland areas, documented ore mining activities in the Tatra Mts (18th and 19th centuries; Landenberger et al. 2003) could have facilitated LDD towards metal enriched sites in Poland (e.g. moving miners or material transport). Such potential anthropogenic seed transport could at least partly explain the admixture signals that we observed within the NM ecotype, as the spontaneous hybridization between M and NM groups seems unlikely in light of the geographic isolation of these two genetically divergent clusters. An alternative explanation for admixed populations is an incomplete lineage sorting, which may still be ongoing due to the longevity of B. laevigata and its capacity to reproduce via vegetative propagation.
As an alternative scenario to recent LDD, ABC results suggest that the divergence between M and NM populations is ancient. Assuming a generation turnover in B. laevigata of 10 to 100 years (Ch. Parisod, personal communication), we roughly dated the divergence time (1,200 generations ago) of M and NM populations between 12 and 120,000 years before present (BP), suggesting that the populations that colonized metalliferous lowland sites may indeed predate the LGM (20-26,000 years BP, Clark et al. 2009). The calcareous slopes and cliffs of the Tatra Mts offer suitable habitats for B. laevigata and probably served as refugia for this and other montane taxa during the glaciations. This is consistent with the interglacial or even pre-glacial relic character of other diploid populations of B. laevigata from Central Europe (Manton 1934;Tremetsberger et al. 2002). The co-occurrence of two genetic clusters of B. laevigata in our study could thus be explained by well-known profound perturbations of the European flora during the last glaciation (Sánchez Goñi et al. 2005), which are considered the principal cause of geographically disjunctive species distributions at both continental and local scales (Bettin et al. 2007;Hewitt 1999Hewitt , 2000Kropf et al. 2006;Taberlet and Cheddadi 2002). This would rather suggest that the occurrence of B. laevigata in the Bolesław -Olkusz region largely predates the development of mining and industrial activities from the Middle Ages. This is in accordance with palaeobotanical studies conducted 50 km from the Bolesław -Olkusz region that suggested that, together with other mountain species, current lowland populations of B. laevigata are remnant of local pre-glaciation populations (Szafer 1930). Under such an assumption, the intriguing disjunctive distribution of B. laevigata in Poland would result from a fragmentation of species range during the Holocene, rather than LDD events.
Biscutella laevigata is regarded as a pioneer and shade-intolerant species poorly competing with other plant species (Dannemann 2000;Oberdorfer 2001). Its ecological affinities have apparently facilitated the colonization of metalliferous habitats, which recently became anthropogenic refugia for lowland populations. Based on our results (genotyping and soil metal concentrations) we assume that contrasted combinations of environmental constraints at metal-polluted sites may result in varying grades of selection pressure and require specific genotypic adaptation. Such site conditions continue to affect establishing populations of B. laevigata. Similar patterns were recently shown in another pseudometallophyte Viola tricolor from the Bolesław -Olkusz area (Słomka et al. 2011), where the existence of two genetic ecotypes among waste heap populations appears to be associated with an unequal horizontal and vertical distribution of metals in the soil.

Conclusions
This study on the genetic structure and divergence of all known M and NM populations of B. laevigata in Poland has shown that the lowland populations are probably relic of previously well-established populations that had a larger distribution range prior to the LGM. This old vicariance between geographically separated lowland (mostly M) and alpine (NM) populations continues to be modified by the ongoing microevolution processes caused by variable selective pressures at mosaic-like metal-contaminated sites. Interestingly, environmental conditions -in particular metal concentrations in the soil -appear to more strongly influence the genetic structure than geographic distance. A significant reduction in the genetic diversity (founder and bottleneck effects) in M compared to NM populations was associated with the colonization of polluted sites and/or evolution of M populations. As a consequence, populations from natural and anthropogenic locations have adapted to different environmental conditions and have genetically diverged. Hence, they should be considered as two distinct edaphic ecotypes. Our study emphasizes the importance of combined natural and anthropogenic historical impacts on the microevolution of M populations, which need to be considered in future investigations of pseudometallophytes. Given the clear patterns of genetic divergence identified in our study, the M and NM populations of B. laevigata from southern Poland represent particularly interesting cases for analyzing the specific contributions of environmental and anthropogenic drivers of microevolution to the emergence of adapted ecotypes.