2.1 Study sites
This study was carried out in the surroundings of the former lead (Pb) and zinc (Zn) smelter named “Metaleurop Nord” located in Northern France (Noyelles Godault, Hauts-de-France, France). The smelter was in activity from 1894 to 2003 and was the only producer of primary Pb in France. Its pyro-metallurgic process had generated large quantities of dust containing many metals. For instance, about 1.0 ton of cadmium (Cd), 16.9 tons of Pb and 31.6 tons of Zn were released in 2002, despite the implementation of technical improvements during the 1970s (DRIRE, 2003). An area of approximately 120 km2 around the smelter is still affected by dust emission released from Metaleurop Nord and from another large Zn smelter named “Umicore” close to Metaleurop Nord (Douay et al., 2008). Several studies have revealed high soil TM contamination in the surrounding soils. For instance, concentrations in agricultural top soils are as high as 21 mg kg-1 of dry soil for Cd, 1132 mg kg-1 for Pb, and 2167 mg kg-1 for Zn (Sterckeman et al., 2002), whereas concentrations of these metals reached up to 67, 4890 and 2685 mg kg-1, respectively, in the upper organic layers of grasslands around the smelter (Sterckeman et al., 2000). In the study of Fritsch et al. (2010), total concentrations of the three metals in soils sampled from wooded patches reached up to 236, 7331 and 7264 mg kg-1 of dry soil, respectively.
The present study was undertaken on six sites of 25 ha (500m x 500m) in the surroundings of Metaleurop Nord. The area is a densely urbanized zone mixed with arable lands (ploughed fields mostly, grassland), forest and planted woods, and shrublands in derelict lands. Five sites among the six ones were located within an area of 40 km2 (5 x 8 km) including the smelter in its center (50°25’42 N and 3°00’55 E). The other site was located about 10 km northeast of the former smelter and is considered as a control site due to its relatively low soil TM contamination: 0.9 - 2.4, 43 - 200, 89-278 mg kg-1 of dry soil for Cd, Pb, and Zn, respectively (Fritsch et al., 2010), which are similar to the regional reference values (Sterckeman et al., 2007). The six sites were chosen for representing a gradient of TM soil concentrations (control, low, moderate, and high contamination, referred to as “TE2”, “103”, “117”, and “097”, respectively, Table 1).
2.2 Field inventory of plants and invertebrates
In each of the six sites, plants, flying and ground-dwelling invertebrates were surveyed as summarized in Figure SI1 and described below, focusing on woody habitats (natural forests, tree plantations such as poplar groves, or copses and hedgerows in natural or cultivated lands and urban parks). Invertebrates were surveyed in woody patches where soils had been sampled and analyzed in previous studies (Fritsch et al., 2010).
Vegetation survey was realized on each of the six study sites from 4th June to 5th September 2012. Vascular plant taxa were identified separately in three different strata: tree stratum (woody species > 8 m high), shrub stratum (woody species < 8 m high), and herbaceous stratum. Plants were identified in the field at species level following Dudman and Richards (1997) and Lambinon et al. (2004). Cover-abundance was visually estimated as vertically projected area for each species following Braun-Blanquet et al. (1952). In total, 236 plant taxa were listed in the three strata: 25 species in tree stratum; 42 species in shrub stratum and 193 species in herbaceous stratum (N.B. some plant taxa were observed in several strata).
Inventory of invertebrates was carried out on each of the six study sites in spring (April) and in autumn (September and October) 2012 by two types of trapping: pitfall traps and yellow pan traps. It is important to note that this sampling aimed at providing a snapshot of invertebrate communities rather than a full description of invertebrate communities, which would have required a longer sampling period.
The pitfall trap is one of the most frequently used methods for sampling epigeic invertebrates such as ground beetles, rove beetles, wandering spiders, and ants (Leather, 2005). Three 800ml polypropylene beakers with neither roof nor preservative fluid were set at 15m intervals between each of them in woody patches where plants were present in all the three strata, constituting and hereinafter referred to as a “trap line”. From two to eight trap lines were set per season (i.e. spring and autumn) on each of the six study sites (i.e. 36 trap lines were used in spring and 32 in autumn, Table 1).
The yellow pan trap is frequently used for sampling flying insects. Trapping color plays a determinant role in the effectiveness with which different insect groups are caught, but yellow color is most efficient for catching a wide range of phytophagous insects and their predators or parasitoids (Kirk, 1984). One to three yellow pan traps with soap mixed water were set on the ground on each of the six study sites and for each season in woody patches where plants were present in all the three strata (i.e. 12 yellow pan traps were used in spring and 12 in autumn, Table 1).
Locations of the two types of traps were geo-referenced. Both pitfall and yellow pan traps were checked every morning for three consecutive days, and then removed from the field. Captured invertebrates were stored in ethanol or in freezer at -20°C and then identified in laboratory at the finest possible taxonomic level by morphological characteristics. The main references used for invertebrate determination were Coulon (2003), Forel and Leplat (2001), Jeannel, (1941) and Trautner and Geigenmueller (1987).
Fauna captured by pitfall and yellow pan traps were considered as “ground-dwelling (GD)” and “flying” invertebrates, respectively. Collembolans were removed from our inventory. Most of the individuals were identified at family level (GD 74% and flying 78%), but other were identified at order (GD 24% and flying 21%) or class levels (GD 2% and flying 1%) (cf. Supporting Information SI Spreadsheet file). Twenty-four taxa at different taxonomic levels (three classes, three orders and 18 families) were listed in the GD invertebrates, while 95 taxa (one class, four order and 90 families) were listed in flying invertebrates. Among the 90 flying invertebrate families, six families of Diptera larvae (Chironomidae, Culicidae, Limoniidae, Psychodidae, Simuliidae, and Sciomyzidae) mainly occur in aquatic habitats, and seven families of Diptera larvae (Ceratopogonidae, Tipulidae, Dolichopodidae, Empididae, Phoridae, Syrphidae, and Scathophagidae) could occur in aquatic or semi-aquatic habitats, depending on genus/species.
2.3 Data preparation
2.3.1 Diversity index choice
Richness (S: number of different taxa) and Simpson’s diversity index (D: 1/∑Pi2, where Pi is the proportional abundance of taxa i) were used to contrasting total number of taxa (richness) to number of abundant taxa at habitat patch level (Jost, 2006). Simpson’s evenness (E: D/S) was calculated and used as another variable with reference to the proportion of dominant taxa among all taxa. Abundance of all taxa (N) was also added as another information about the community. The four indices are hereinafter referred to as “alpha diversity”.
Spatial variation in composition among communities (i.e. beta diversity) was estimated by using the total variance of the site-by-taxa community data (Legendre et al., 2005; Legendre and De Cáceres, 2013). The beta diversity can be partitioned into two matrices representing “replacement” and “richness difference” (Borcard et al., 2018), and each matrix can be analyzed in relation to explanatory environmental variables (Legendre and De Cáceres, 2013). For all types of plant strata and of invertebrates, dissimilarity matrices for beta diversity were built from presence-absence of each taxon in communities because binary dissimilarity coefficients produce more relevant results than quantitative indices when taxa are largely different among communities (Legendre, 2014). The two matrices for replacement and richness difference were hereinafter referred to as “beta diversity”.
2.3.2 Calculation of diversity indices and matrices.
Alpha and beta diversity of invertebrates were based upon individuals captured by 68 trap lines for GD invertebrates and 24 yellow pan traps for flying invertebrates. Diversity of plants were measured for each stratum (tree, shrub and herb) based on cover-abundance (m2) of each species present in an area of 1000 m2 around trap lines. As pitfalls and yellow pan traps were not precisely set at the same locations, alpha and beta diversity of plants were measured for each type of invertebrate traps. Those areas are hereinafter referred to as “buffers”. As the plant inventory was done once between June and September, their presence and relative cover-abundance were considered to be similar at the two seasons for further statistical analysis.
Both alpha diversity and beta diversity were calculated at species level for plants. Calculation of diversity was carried out at family level for invertebrates when available, otherwise at order or class levels (cf. SI Spreadsheet file). Dissimilarity matrix for beta diversity was built using the Jaccard dissimilarity coefficient and partitioned into matrices for replacement and richness difference, hereinafter referred to as “replacement” and “richness difference”, respectively.
2.3.3 Soil data
Soil properties and soil TM concentrations were referred to Fritsch et al. (2010). As concentrations of Cd, Pb and Zn in soils were highly correlated (Pearson’s r > 0.9, p-value of correlation test < 0.001), only Pb concentration in soil was used as a proxy of soil TM contamination in our statistical analyses. Soil pH and organic carbon content (g kg-1) in soil, considered as a proxy of the organic matter (OM) content in soils, were used as soil properties importantly related to metal bioavailability (Bradham et al., 2006; Giller et al., 1998; Visioli et al., 2013). Soil pH was also positively correlated with soil TM concentrations (Spearman’s rho > 0.6, p-value of correlation test < 0.001 for the three TM). For each buffer, we used soil TM contamination and soil properties of the nearest soil sampling point from the given buffer (i.e. no more than 50 m). Trace metal soil contamination (as represented by Pb soil concentrations), OM content and pH were hereinafter referred to as “Pb”, “pH”, and “OM”, respectively, and linked to each buffer.
2.4 Statistical Analysis
2.4.1 Data transformation
Lead concentrations and OM, as well as abundance (i.e. total cover-abundance of plants or total number of individual invertebrates), were logarithmically transformed because of their skewed distributions. Alpha diversity and the soil variables were then scaled to zero mean and to unit variance for each variable because of their different unit.
2.4.2 Relationships between plant diversity and soil properties
Before assessing the effects of plant diversity on invertebrate communities, the relationship between plant diversity, soil TM contamination and soil properties was explored. Both alpha diversity indices and beta diversity matrices for plants were calculated in vegetation patches near soil sampling points of Fritsch et al. (2010). The analysis used vegetation around 17 soil sampling points, where plants were available for all three strata.
A redundancy analysis (RDA) was executed for the alpha diversity in relation to soil TM contamination and soil properties. A forward selection of significant explanatory variables was carried out (Borcard et al., 2018). Proportion of variance explained by the selected explanatory variables was indicated by an adjusted R2 (R2adj) (Peres-Neto et al., 2006). Relationships between plant beta diversity, soil TM contamination and soil properties were assessed by using the distance-based RDA (dbRDA) (Legendre, 2014). Briefly, a principal coordinate analysis (PCoA) was carried out for each dissimilarity matrix after square-rooted transformation. Their principal coordinates were used as response variables and a forward selection of significant explanatory variables was carried out, and R2adj was measured. Furthermore, species presence-absence data were a posteriori projected on the ordination plot using weighted averages (Borcard et al., 2018), which shows how occurrence of species in communities is affected by the environmental factors analyzed.
2.4.3 Relationships between diversity of invertebrates, diversity of plants and soil properties
The partial RDA (pRDA) and the variation partitioning was applied for invertebrate alpha diversity, using soil TM contamination, soil properties and plant alpha diversity as explanatory variables (Borcard et al., 2018). The pRDA for plant diversity was executed as follows: after RDA for invertebrate alpha diversity data in relation to soil TM contamination and soil properties, the residual variation of this RDA (i.e. variation of invertebrate diversity data non-explained by soil properties) was handled by another RDA in relation with plant alpha diversity. This was vice versa for the pRDA in relation to soil properties. Variation explained by selected variables of each explanatory matrix, as well as variation explained jointly by them, were shown using Venn diagram. If one of the two explanatory matrices was not significantly related to the response matrix, ordinal RDA was carried out. The partial dbRDA was applied for invertebrate beta diversity, using soil TM contamination and soil properties and plant alpha diversity as explanatory variables. Invertebrate presence-absence data were then a posteriori projected on the ordination obtained.
All statistical analyses and graphics were performed by the statistical software R (ver. 3.6.1; R Development Core Team). PCA, PCoA and RDA were handled with the “vegan” package. The function “forward.sel” of the package “adespatial” was used for forward selection. The functions “beta.div.comp” and “dbRDA.D” from Legendre (2014) were used for building replacement and richness difference dissimilarity matrices and for carrying out accurate significance test for dbRDA, respectively.