Characterization of the Prayon site
The two hyperaccumulators A. halleri and N. caerulescens colonize the Prayon site (near Liège, Belgium, Fig. S1), enabling the comparison of their hyperaccumulation capacities in a natural location and the study of the specificities of the microbial composition of their rhizospheres. Shoot, rhizosphere, bulk 1 and bulk 2 (top and bottom of the core) samples of each plant species were collected in six plots on this calamine site (Fig. 1) in order to determine if their convergent evolution was reflected in the rhizosphere microbiota.
pH, water content, temperature and element composition, including zinc and cadmium, were measured in bulk 1 and bulk 2 soil samples (Fig. 2). pH and temperature were stable through the six plots: temperature ranged between 17.8 and 19.2°C and pH between 6.04 and 6.59. This pH is lower than in other studies of A. halleri (Corso et al., 2018; Stein et al., 2017) and some of N. caerulescens (Luo et al., 2018; Rees et al., 2019; Rosenfeld et al., 2018) even if it seems more common for N. caerulescens (Luo et al., 2018; Rosenfeld et al., 2018). Water content was slightly more variable, ranging between 29.52 and 40.47% of fresh weight for bulk 1 and between 27.41 and 58.14% for bulk 2, which is coherent with previous study (Luo et al., 2018).
Overall, a limited variation is observed for extractable (HCl) and exchangeable (BaCl2) element concentrations between soil cores associated with A. halleri or N. caerulescens. Between bulk 1 and bulk 2 samples, only a few extractable elements varied (Al, Fe and Na), in contrast to exchangeable elements. Fe and Na are the most variable elements across plots (Fig. S4). As expected based on the history of the site (Denaeyer-De Smet and Duvigneaud, 1974; Duvigneaud and Jortay, 1987; Ramaut, 1964), a large zinc pollution is measured. Current values are largely higher than those measured in 197433, even with respect to exchangeable rather than extractable zinc. This could be explained by the differences in the techniques used since the activity of the factory stopped in 1970.
Moreover, extractable Zn concentrations (i.e., 4,934 to 25,051 mg kg-1 for bulk 1 and 6,326 to 13,829 mg kg-1 for bulk 2) are largely superior to the reference values for Wallonia soil (“Législation/décret gestion et assainissement sols,” n.d.). Indeed, natural soil is expected to contain a maximum Zn concentration of 196 mg kg-1 (DW) and the threshold for necessary intervention on industrial soil set by Wallonia regulations is 3,000 mg kg-1. The same is true for cadmium, for which the intervention threshold is fixed at 20 mg kg-1. Here in Prayon, even the exchangeable concentrations (42.45 - 174.77 mg kg-1) are largely above this limitation. This is not the case for Pb and Cu, which have high values (extractable concentrations are respectively 500-1500 and 250-625 mg kg-1) but remain below intervention thresholds (respectively 1840 and 600 mg kg-1).
If we compare with concentrations in natural soils reported in other studies, Zn, Cu and Cd exchangeable concentrations are in the higher range of contaminated sites in which A. halleri can be found, and Pb exchangeable concentration is also within metallicolous site ranges (Stein et al., 2017). Moreover, Zn and Cd extractable concentrations are above most natural metallicolous sites colonized by A. halleri. Concerning N. caerulescens, the Prayon values of exchangeable Cd and Zn surpass total soil concentration found in other studies of natural soil (Rees et al., 2019) except one that is in the same range (Rosenfeld et al., 2018). This could explain the low plant diversity found in the sampling site. Along with A. halleri and N. caerulescens, only Armeria maritima, Viola lutea and Festuca ovina were identified, three metal hypertolerant but not hyperaccumulator plants.
Hyperaccumulation in plant shoot
Soil parameters influence the accumulation levels in shoot, thus shoot samples (composed of five pooled individuals) were collected within six 2-m2 plots for each species (Fig. 1) and submitted to ionome profiling and quantification of metal transporter gene expression.
Micro- and macronutrient shoot concentrations were mostly similar between the two plant species, with the exception of Mg, P and Mn concentrations, which were higher in A. halleri than in N. caerulescens (Fig. 3). Metal accumulation in shoot differed significantly between the two species only for Cd, where N. caerulescens hyperaccumulated more. In contrast, Zn and Pb concentration ranges did overlap between shoots of both plants. Finally, A. halleri and N. caerulescens hyperaccumulated both Zn and Cd above the defined hyperaccumulation thresholds of 3,000 and 100 mg kg-1 (DW), respectively (Krämer, 2010). Even if N. caerulescens is known to hyperaccumulate Cd, it is not the case for all ecotypes. It is particularly the case of the Prayon ecotype, which was described as unable to hyperaccumulate Cd (Rosenfeld et al., 2018). However, most studies used soils with lower Cd concentration, therefore limiting Cd accumulation. Furthermore, shoot Cd accumulation in the Prayon site is much lower than observed in the Ganges ecotype of N. caerulescens, a bona-fide Cd hyperaccumulator ecotype (Halimaa et al., 2019).
We next examined the expression levels of metal transporter genes that were previously shown to be highly expressed in the shoot of A. halleri and N. caerulescens and linked to hyperaccumulation and/or tolerance (Hanikenne et al., 2008; Talke, 2006; van de Mortel et al., 2006): HMA3 (Zn and Cd vacuolar storage (Talke, 2006; Ueno et al., 2011)), HMA4 (Zn and Cd xylem (un)loading and distribution (Craciun et al., 2012; Hanikenne et al., 2008)), ZIP4/ZNT1 (Zn uptake and distribution (Lin et al., 2016; Talke, 2006)), ZIP6 (Cd tolerance (Spielmann et al., 2020; Wu et al., 2009)) and MTP1 (Zn vacuolar storage (Dräger et al., 2004; Gustin et al., 2009; Shahzad et al., 2010)). Note that information on the expression of those genes in plants in the field is very limited. The HMA3, HMA4 and ZIP4/ZNT1 genes showed the largest variation in expression level (up to 3.4-fold) throughout samples, both in A. halleri and N. caerulescens. MTP1 seemed also to vary but to a lesser extent (Fig. S5, S6). Correlations of gene expression were computed with the element concentrations in shoot tissues. In A. halleri, a relation between the Pb concentration and ZIP4 expression is observed (p-value=0.0069), as well as between P and MTP1 (p-value=0.047). In N. caerulescens, a correlation is observed between the Mn concentration and the HMA4 expression (p-value=0.00046) and between Mg, Pb and K and MTP1 (p-value=0.02, 0.035 and 0.062, respectively), as well as Pb and HMA4 (p-value=0.052) (Fig. S7).
Correlations were also observed between the soil concentration of some elements and the expression of related transporter genes. In A. halleri, it is the case between the concentration of Zn and the expression of ZIP4 (p-value=0.0099) as well as between Cd and HMA4, Cd and ZIP6, Na and ZIP4 (p-value=0.019, 0.031 and 0.049, respectively). In N. caerulescens, correlations were found between ZNT1, HMA3, ZIP6 and HMA4 and several elements (Fig. S8).
Those correlations are relatively weak and at this point, disentangling causal relationships between gene expression and metal concentrations would be hazardous. Differences in the observations in A. halleri and N. caerulescens might be explained by their different physiology (Merlot et al., 2018), including distinct metal storage sites, in mesophyll and in epidermis for A. halleri and N. caerulescens, respectively.
Interestingly, in spite of a generally high level of zinc and cadmium pollution, the heterogeneity of the site of Prayon seems to correlate with the expression of some metal transporter genes. However, no correlation was found between the amount of zinc present in shoot and the amount present in soil, in contrast to Stein et al. 2017 (Stein et al., 2017) (Fig. S9).
On the optimal sequencing depth
Since microorganisms present in the rhizosphere impact plant physiology and can be species-specific, a major goal of this study was to characterize as accurately as possible the microbial diversity of our soil samples, i.e., the number and relative abundance of each OTU or taxon. To this end, DNA was extracted from all soil samples (rhizosphere and bulks) and SSU rRNA (16S/18S) gene regions were amplified either with universal primers or with bacterial and archaeal-specific primer pairs for corroboration. Amplicons were sequenced using the Illumina MiSeq technology.
We first tested whether the planned sequencing depth enabled the identification of the whole microbial diversity in our samples. This pilot analysis was conducted on three samples (rhizosphere, bulk 1 and bulk 2) from five pooled A. halleri individuals sampled in the same plot. Five different depths were tested, obtained by in silico sub-sampling of the available reads, resulting in a relative depth from 1x to 12x in each of the three samples. Reads from each depth were treated independently and processed through the pipeline described in Material and methods to identify the OTUs (i.e., potential organisms) present.
As observed in Fig. S10, the number of different taxa discovered for any sequencing depth and any taxonomic level reached a different plateau. Hence, at deeper depth, more organisms are found when replaying the whole pipeline, suggesting that these rarefaction curves are not informative about our ability to discover (or not) all the microbial richness in the samples. Because we observe this phenomenon whatever the taxonomic level analyzed, it is unlikely to be (entirely) due to “new” sequences resulting from undetected chimeras or sequencing errors. Moreover, similar results were obtained with all primer pairs.
Yet, as many OTUs (up to 12.7% at the deepest depth) actually correspond to unclassified organisms according to SILVA taxonomy, it was difficult to determine whether increasing the sequencing depth really allowed sampling a significantly larger biodiversity. To address this issue, we computed a phylogenetic tree including all the discovered OTUs and labelled the leaves with the depth(s) at which each OTU was recovered (Suppl. File at [https://doi.org/10.6084/m9.figshare.14376785.v1]). This tree was then used to build a pseudo-taxonomy with five nested levels of pseudo-taxa accounting for all OTUs, including those unclassified in the SILVA taxonomy. These analyses revealed that new pseudo-taxa consistently appear at all levels when increasing depth (Fig. S11). Hence, between 16.4 and 22.7% of the clusters segmenting the phylogenetic tree are only observed at the two deepest sequencing depths.
Those results suggest that there is no good way to ascertain the optimal sequencing depth to uncover the whole microbial richness, even if some studies try different approaches to resolve this problem (Alteio et al., 2020; Fuks et al., 2018; Xiao et al., 2018). Therefore, comparison between studies or samples should take into account the depth parameter in addition to the primer pair, the sequencing method and the bioinformatics pipeline used. In the present work, the sequencing depth was eventually set at 4x the depth of the pilot study, as it was the minimal depth providing a reasonable account of the real diversity.
Assembly and quantification of SSU rRNA OTUs
Altogether, a total of 36 samples (rhizosphere and bulks from both species in the six locations), along with extraction and amplification blanks, were sequenced by Illumina after amplification of a SSU rRNA fragment. The results presented in the following were obtained with the universal primer pair, which uncovered the largest richness (compared to bacterial and archaeal primer pairs). Indeed, data obtained with bacterial and archaeal primer pairs had to be normalized at 1,000 reads (instead of 30,000; see below). Moreover, they were not especially Bacteria or Archaea-specific, contrary to our expectations, and did not allow us to find additional organisms beyond those obtained with the universal primer pair. The corresponding data is nevertheless available at [https://doi.org/10.6084/m9.figshare.14258528.v1].
Reads were analyzed through the bioinformatics pipeline described in Material and methods in order to generate a normalized OTU table in which each OTU has been assigned a taxon based on SILVA taxonomy. The normalization was fixed for all samples at 30,000 reads. It was the minimal number of reads for which we could have a reasonable taxonomic diversity with a maximum of our samples. The resulting OTU table counts 9,796 OTUs (out of the 11,328 originally identified by usearch) and 30 out of the 36 samples were conserved. Because of the 30,000-read normalization, we have lost the bulk 2 sample from A. halleri in the second plot and five samples from the sixth plot (i.e., all except the bulk 1 from N. caerulescens). Finally, after the filtration step, 9,323 OTUs remained and were used in the following analyses.
Microbial diversity at the OTU level
The taxonomic diversity was then compared between these 30 samples. A non-metric multidimensional scaling (NMDS) was first computed to visualize whether samples had similar (or different) microbial diversity comprising the number of different organisms and their abundance (Fig. 4a). A clear separation between each part of the core is observed, i.e., between rhizosphere, bulk 1 and bulk 2. It also shows a gradient of diversity from the roots to the bottom of the core, i.e., rhizosphere then bulk 1 (top of the core) then bulk 2 (bottom of the core). Interestingly, the projected distance between rhizosphere and bulks is larger than between the two bulks, suggesting a more different microbial diversity in the rhizosphere. Moreover, the variance within each type of sample apparently follows the same gradient.
Since the major variability, i.e., differences found in the diversity between samples, seems to be driven by the physical distance from the roots, which is in line with other studies (Praeg et al., 2019; Song et al., 2020), this analysis could hide potential differences in the microbial diversity between the rhizospheres of the two species (A. halleri and N. caerulescens). To check that possibility, another NMDS was computed for those specific samples (Fig. 4b). While samples from either plant visually form distinct groups, their separation is not very clear, suggesting only subtle differences between rhizosphere compositions, less important than the differences between sample types (i.e., rhizosphere, bulk 1 and bulk 2). This is coherent with the distribution of OTUs in Venn diagrams, according to sample type or plant species (Fig. S12).
Statistical values for alpha- and beta-diversity were then computed at the OTU level within and between samples, respectively (Fig. 5) to better understand the variations revealed in the NMDS plots. A decreasing gradient of richness is observed from rhizosphere to bulk 2 (Fig. 5a), i.e., the number of different organisms is the highest in rhizospheres, which is fully coherent with the first NMDS analysis (Fig. 4a).
However, Simpson’s index of alpha-diversity, which takes into account both the number of “species” present and their abundances, highlights that the overall variability in bulk 1 is similar to the one found in the rhizosphere (Fig. S13a), even if there are more organisms that are different in the rhizosphere (Fig. 5a). Simpson’s index also shows a larger variability between bulk 2 samples (Fig. S13a), which is confirmed by Euclidean distances of beta-diversity, taking abundances into account too (Fig. 5b). Similarly, beta-diversity measures also show that bulk 2 differs more from bulk 1 and rhizosphere, which are very similar in terms of Euclidean distances computed on OTUs.
In summary, the taxonomic diversity is mostly driven by the physical distance from the roots, along a decreasing gradient, rhizospheres having the highest richness and Simpson’s diversity but a smaller variability between samples than bulks, with bulk 2 having the highest variability and the smallest number of different organisms. Thus, the rhizosphere appears to be a more favorable microenvironnement than bulks but with more constraints. However, those constraints appear similar between the two plant species, considering the convergence of the microbial communities observed in their rhizosphere.
Microbial diversity at different taxonomic levels
The microbial composition (top-10 taxa) was then compared between samples at different taxonomic levels (Fig. 6a). Differences appear between rhizosphere and bulks, and maybe between plant species too. At the phylum level, even if the most represented phylum is Proteobacteria for all samples, it is less abundant within bulk 2. More Bacteroidetes and Actinobacteria are present in the rhizosphere. Moreover, Acidobacteria, Chloroflexi and Firmicutes are less represented in the rhizosphere than in bulks, which may depend on metal contamination levels as suggested by observed results in the presence of severe Cd-pollution in a previous study (Song et al., 2020). Interestingly, unclassified taxa (according to SILVA taxonomy) are among the 10 most represented taxa at each taxonomic level.
Based on the results obtained with OTUs (Fig. 5a), we expected to observe a larger richness in the rhizosphere. However, if we take into account the “other” taxa that are clustered together (i.e., the remaining taxa out of the 10 most abundant), it does not seem to be the case at the phylum level (Fig 6a). Indeed, bulk 2 has not only more “other” phyla but also a larger proportion of unclassified and unknown (= undef_tax) organisms. This difference might be biologically genuine and explained by visualization at a higher taxonomic level (i.e., phylum) than the OTUs used previously. As a case in point, the richness shows a loss of signal and a large variability that is not in line with observations on OTUs throughout the different taxonomic levels (Fig. S14). Regarding beta-diversity measures, they are even more variable across levels (Fig. S15).
Alternatively, the systematically large proportion of unclassified and unknown organisms might also impact the analyses, considering that information loss is greater at lower taxonomic levels, i.e., it is more difficult to reliably classify an OTU at the genus level than at the class or phylum level (Table S4, Fig. 6a). In detail, 6,186 out of 9,323 identified OTUs (66%) are classified at any level, and only 1,794 (19%) are affiliated down to the genus level. The remaining OTUs (15%) are classified at various higher taxonomic levels, including 350 OTUs (4%) above the phylum level, thus remaining unknown for most practical purposes. Furthermore, taxonomic affiliation levels are heterogeneous between sample types, with more reads classified at the genus level in the rhizosphere than in bulk 2 and the opposite at the phylum level (Fig. S16, Fig. 6a).
Microbial diversity with pseudo-taxonomic affiliations
To bypass these biases created by missing taxonomic affiliations, a pseudo-taxonomy was computed following the same strategy as for the sequencing depth pilot study (see above), but this time using all OTUs remaining after filtration (9,323, Fig. S2) to compute the phylogenetic tree (Table S4). As shown in Table S5, the number of different organisms found in each pseudo-taxonomic level is constant.
With this pseudo-taxonomy, richness measures were computed anew (Fig. S17). No particular signal seems to emerge at pseudo-class and -phylum levels. However, the variability pattern observed at pseudo-species, -genus and -family levels is congruent with the pattern found with OTUs (Fig. 5a). Indeed, the richness is higher in the rhizosphere and follows a gradient of decreasing diversity through sample types. This was not observable above the OTU level when using the SILVA taxonomy (Fig. S14).
Beta-diversities were then recomputed on this pseudo-taxonomy (Fig. S18). They appear coherent both with Simpson’s index (Fig. S19) and earlier analyses carried on OTUs (Fig. 5b). Indeed, there is a gradient throughout pairwise comparisons from rhizosphere to bulk 2 at the pseudo-order, -class and -phylum level, thereby confirming a smaller variability in the diversity found in the rhizosphere than in bulks, with bulk 2 being the most variable. Moreover, beta-diversity also shows a significant difference between the rhizospheres from A. halleri and N. caerulescens (Fig. S18) at multiple levels, as suggested by the second NMDS analysis (Fig. 4b) and the composition found with the pseudo-taxonomy (Fig. 6b). Indeed, differences (Fig. S18) between rhizosphere composition of the two plant species at pseudo-phylum, -class and -order are visually observed and some variability appears at pseudo-family and -species levels.
Euclidean distances between rhizospheres of A. halleri are almost systematically lower than those between rhizospheres of N. caerulescens (Fig. S18). This may suggest a sharper selection caused by the roots of A. halleri than by those of N. caerulescens.
Differences in the microbial diversity within rhizospheres
To investigate a possible difference of microbial composition in rhizosphere samples from both plant species, a paired t-test was computed across all six plots for all pseudo-taxonomic levels (Table S6a), using only the most abundant pseudo-taxa (i.e., those shown in Fig. 6b). However, not a single pseudo-taxon abundance appeared as significantly different at any taxonomic level, the best candidates being cluster 41 at the pseudo-family level and cluster 536 at the pseudo-species level, with a corrected p-value of 0.08 and 0.06, respectively. Cluster 41 is composed of OTUs that were taxonomically associated to Chloroflexi (76%) and Armatimonadetes (4%) taxa (both from Terrabacteria group), as well as unclassified OTUs (20%), whereas cluster 536 is composed of OTUs assigned to Pirellulaceae (96%) and of a few unassigned OTUs (4%). Interestingly, those phyla are often cited when soil conditions vary (Edwards et al., 2015; Mukhtar et al., 2018; Sánchez-Marañón et al., 2017; Song et al., 2020; Sun et al., 2018; Wu et al., 2019), except Pirellulaceae, which might be an indicator of soil pH (Hermans et al., 2017).
In contrast, a comparison between bulk 1 and rhizosphere of A. halleri identifies 2-5 pseudo-taxa with a significant difference at each pseudo-taxonomic level (Table S6b), as does the same comparison for N. caerulescens samples (Table S6c). Among all those differences, Gemmatimonadetes, also varying in other studies (Sun et al., 2018; Yue et al., 2020; Zhao et al., 2020), Chloroflexi and Armatimonadetes (along with unknown and undefined organisms) are found in several clusters at various levels for A. halleri (cluster 435 at the pseudo-genus level, clusters 71 and 41 at pseudo-family, cluster 1 at pseudo-order and cluster 1 at pseudo-class). Cluster 536 at the pseudo-species level highlights again a difference in Pirellulaceae. Concerning the most significant differences found between bulk 1 and the rhizosphere of N. caerulescens, Chitinophagaceae compose the whole cluster 418 at the pseudo-species level. This phylum was previously identified as an indicator of aluminium concentration in soil (Hermans et al., 2017), but it was also shown to be mostly found in rhizosphere (Song et al., 2020).
These statistical analyses are globally in line with the results obtained at the OTU level (Fig. 4, Fig. 5) and corroborate the insights brought by our alpha- and beta-diversity analyses based on pseudo-taxonomy (Fig. S19, Fig. S18). While they demonstrate that relying on a pseudo-taxonomy can indeed be enlightening in a study like ours, they also confirm that the microbial communities of the two plant species studied in this work are essentially convergent. Such an outcome might be explained either by the spatial proximity of the roots sampled in the field and/or by the phylogenetic proximity of the corresponding plant species. Indeed, these two species belong to Brassicaceae and their rhizospheres might exert a similar selection on a possibly similar pre-existing microbial community. A sampling of rhizospheres from non-hyperaccumulator plant species present on the same site could be useful to ascertain this hypothesis, such as Armeria maritima or Viola lutea. Another explanation could stand for the necessity of taking into account less abundant taxa. This would be in line with statistical differences obtained in beta-diversity (Fig. S18) and not with major taxa. Moreover, less abundant unknown and undetected taxa may be of importance in the specificity of microbiomes (Thomas and Segata, 2019).