Metal-induced bacterial interactions promote diversity in river-sediment microbiomes

Background Anthropogenic metal contamination results in long-term environmental selective pressure with unclear impacts on bacterial communities, key players in ecosystem functioning. Since metal contamination poses serious toxicity and bioaccumulation issues, assessing their impact on environmental microbiomes is important to respond to current environmental and health issues. Despite elevated metal concentrations, river sedimentary microbiome near the MetalEurop foundry (France) show unexpected higher diversity compared to upstream control site. In this work, a follow-up of the microbial community assembly during a metal contamination event was performed in microcosms with periodic renewal of the supernatant river water. Results Sediments of the control site were gradually exposed to a mixture of metals (Cd, Cu, Pb and Zn), in order to reach similar concentrations to MetalEurop sediments. Illumina sequencing analysis of 16S rRNA gene amplicons was performed. Metal resistance genes czc A and pbr A, as well as IncP plasmid content were assessed by quantitative PCR. The outcomes of this study support previous in-situ observations showing that metals act as community assembly managers, increasing diversity. Conclusion This work revealed progressive adaptation of the sediment microbiome through the selection of different metal-resistance mechanisms and cross-species interactions involving public good providing bacteria co-occurring with the rest of the community. pbrA, sequences Cupriavidus metallidurans CH34, Klebsiella Enterobacter cloacae using the ClustalX. consensus used templates for the design of degenerated primers. physicochemical properties of each pair of primers compared Scientific Multiple Primer Analyzer tool [80]. best candidates S5). specificity of primers tested using BLAST (NCBI) and classical PCR of pure Cupriavidus metallidurans CH34 and Pseudomonas putida KT2440 cultures. For quantitative PCR (Q-PCR), the DNA samples were purified using a QIAquick® PCR Purification (QIAGEN) kit to eliminate potential PCR inhibitors. A StepOnePlus Real-Time PCR System (Thermo Fisher Scientific) was used for the quantification of oriT in IncP, czc A and pbr A. The relative concentrations of these genes were calculated as the ratio of the expression of the target gene and that of the reference 16S rRNA gene using the universal primers 518R (5’-ATTACCGCGGCTGCTGG–3’) and 341F (5’- CCTACGGGAGGCAGCAG –3’). The results are presented as a double ratio of the expression of the target gene (IncP oriT, czc A or pbr A) and the reference gene (16S rRNA) at a given time in comparison with the expression at T0. The ratios were calculated according to equation (1) [81]

The use of metals in anthropogenic activities increased their dissemination in the environment [11,12], including in river ecosystems [13]. While some metals are essential nutrients, all have harmful toxic effects at high concentrations. Among prokaryotes, the most widespread metal resistance systems found in in-situ bacterial communities are in the periplasm or anchored in the cytoplasmic membrane [14], and are largely represented by efflux pumps. In addition to the direct effects of metal contamination on biological systems, metal-induced stresses co-select for antibiotic resistance mechanisms, for example by co-selection of mechanisms involved in efflux systems [15][16][17][18].
Furthermore, public good providing bacteria were shown to provide interspecific benefits in a complex community [19]. These bacteria, also called facilitators, detoxify sediments (e.g. via metal precipitation and sequestration) and provide safe micro-niches for metal-sensitive bacteria (e.g. biofilm formation) [20]. Metal contamination was found to maintain [1,21], decrease [22] or increase microbiome richness and evenness [7]. These contrasting observations might be linked to both the site properties (e.g., initial composition of the microbiome, the nature of contaminants, exposure times or local geo-physicochemical parameters).
The Deûle River (Northern France) flows next to the former foundry MetalEurop where sediments are up to 30-fold more contaminated with metals than an upstream control site (the Sensée canal, Férin) [1]. Interestingly, MetalEurop sediments exhibited a highly diverse microbiome compared to the control site [1,7]. Jacquiod and colleagues suggested that metals may play a key role in the assembly of in-situ microbial communities through the coalescence (i.e. the gathering of microbiomes from different ecosystems [23]) of indigenous bacteria and bacteria from different upstream environmental sources (e.g. WWTPs or agricultural sources). Authors suggested that metals increased diversity by (i) preventing the development of metal-sensitive opportunists and (ii) selecting metal-resistant bacteria providing public-goods for the rest of the community. In addition, Horizontal Gene Transfer (HGT) has been suggested to be involved in metal-resistance gene dispersal in the community [7]. If plasmids present in the environment, especially IncP plasmids, have been shown to carry metal resistance genes [24,25], the impact of metals on their persistence and dispersion in the environment is not clear.
In the present study, we incubated control sediments in microcosm for more than 6 months, some of the microcosms being exposed to metals that gradually reached those found in nearby contaminated site MetalEurop: Cd (~ 40 mg/kg), Cu (~100 mg/kg), Pb (~900 mg/kg) and Zn (~3000mg/kg) [1,7].
The overlaying water was renewed periodically with fresh river water to induce community coalescence phenomenon such as observed at natural conditions [6,7]. We performed a follow-up of the taxonomic and functional evolution of river sedimentary microbiomes from the Sensée canal (Férin, France) in the presence of metal concentrations. The community assemblage in sediment was monitored over time using 16S rRNA gene amplicon sequencing and quantitative PCR to follow the metal resistance genes czcA and pbrA, encoding metal efflux pumps, but also, the evolution of the IncP plasmid family.

Results
The α-diversity of river-like metal-impacted microcosms remains high over time.
Species richness and Shannon's diversity index decreased slowly during the 6.5 months of incubation ( Figure 1). Both indices were significantly higher in the control microcosms than in the metal-treated microcosms after 15 days of incubation. Richness fluctuated at the same level in all microcosms, while after 3.5 months, Shannon index was significantly higher in the metal-treated microcosms than in the controls ( Figure 1). Community succession in control and metal-impacted microcosms PERMANOVA revealed significant changes in ß-diversity over time (Table S3). Time, metals and the interaction of both significantly affected the sediment microbial structure. Phylogenetic indices changed over time in both control and metal-treated microcosms: (i) At the start of the experiment, NRI index of the metal-treated microcosms increased and then fluctuated around 0. In control microcosms NRI index immediately decreased and reached a value of -2.20 after 4.5 months but came back to initial values after 6.5 months ( Figure 2A); (ii) β-NRI fluctuated around 0 in the metaltreated microcosms and got close to 2 in control microcosms before coming back to initial values ( Figure 2B); (iii) NTI was high (>2) and relatively stable over time in both conditions ( Figure 2C) and; (iv) β-NTI index started at -1 at the beginning of the experiment and progressively increased to over 0 after 3.5 months ( Figure 2D). Before the first 3.5 months, there were no significant differences in phyla abundance between control and metal-treated microcosms (Table S4). From 3.5 months, significant differences were observed for several phyla, including α-Proteobacteria, β-Proteobacteria, Cyanobacteria and δ-Proteobacteria, Firmicutes, Euryarchaeota, Chloroflexi or Saccharibacteria. In the control microcosms, Cyanobacteria, Euryarchaeota and Saccharibacteria were found to be enriched.
Dynamic interaction network reveals links between keystone species and facilitator bacteria.
The eLSA approach was used to identify simultaneous and delayed OTU occurrences within the microbiome. The clustering coefficients (i.e. degree to which nodes grouped together) were 0.23 and 0.28 and the average number of neighbors was 2.9 and 8.3 for the control and the metal-treated microcosms respectively. In control microcosms, 3 keystone OTUs (i.e. OTU 356 (Hyphomicrobiaceae),, OTU 308 (Chromatiales) and OTU 358 (Verrucomicrobia from Subdivision 3))) were found to be highly connected with the rest of the community ( Figure S4A). In metal-treated  Figure 5). The pbrA level was always significantly higher in the metal-treated microcosms than in the control microcosms ( Figure 5).

Discussion
The community diversity in metal-impacted river-like microcosms Metals affected the community structure after 15 days of exposure by causing a significant decrease of alpha-diversity indices. However, after 3.5 months, this trend was reversed, confirming previous insitu observations [7]. Interestingly, in a preliminary similar experiment, where sterilized water renewal was used, Shannon index in metal-contaminated microcosms remained lower than in control condition ( Figure S5). In the present research, the evolution of phylogenetic indices over time was assessed assuming that: (i) overdispersion (decreased NRI or NTI [26]) would result from diversification of the community in line with the community coalescence process (i.e.,, the constant arrival of new bacteria by water renewal) and (ii) increased relatedness (increased NRI or NTI [26]) would result from environmental filtering (i.e.,, killing or inhibition of some bacteria by metals). In control microcosms, the NRI index decreased under the threshold value (-2), highlighting a phylogenetic overdispersion of the community. This may be explained by a community coalescence process with novel bacteria introduced during water renewal. In metal-treated microcosms and at the very start of the experiment, NRI was found to increase, showing initial environmental filtering most likely for metal-tolerant bacteria, but stabilized around 0 in the following months. Despite metals inhibited overdispersion, the coalescence process may indeed avoid strong environmental filtering. In both control and metal-treated microcosms, β-NTI was first negative but increased over 0 after 3.5 months showing no metal impact in the selection of the nearest phylogenetic neighbors. These results highlighted that the periodic arrival of new bacteria in a metal-contaminated environment helped the microbiome to persist, leading to a fast resilience and an unexpectedly high diversity.

Mechanisms of resilience in metal-impacted microbiome
Previous studies suggested (i) broad-host-range plasmids [7,27] [1,7], and (ii) bacterial cooperation via facilitator bacteria [7], to play an important role in the resilience observed in long-term stressed communities, such as those in metal-impacted river sediments. Broad host range plasmids of the IncP family have been found in areas contaminated with manure, wastewater of industrial origin or river sediments and often carry catabolic, antibiotic-and metal-resistance genes [28,29]. In addition, IncP plasmids were shown to be enriched in metal contaminated MetalEurop sediments, suggesting their role in the long-term resilience of the community [27]. In the present study, IncP plasmid copy depletion in metal-treated microcosms may be explained by the cumulative costs of metal-stress metabolic response, plasmid maintenance, expression and conjugation [27].
Evidences of facilitator bacteria were found in groups of OTUs that responded similarly to the experimental conditions over time (TRGs). Metals were shown to prevent the development of opportunistic anaerobic, acetolactic methanogens Euryarchaeota (mainly the genus Methanotrix),,known for their metal sensitivity [30] (TRG A). Such growth limitations were previously observed in-situ [7]. In addition, metals progressively favored the development of metal-resistant facilitator bacteria that precipitate metals or form large biofilms in the sedimentary structure (B TRGs). For example, in TRG B1 and 3, members of Zoogloea and Dechloromonas were previously found to be highly represented in copper-contaminated activated sludges [31]. Zoogloea has been reported for its propensity to produce massive EPS matrices [32] and to accumulate metal ions [33].
Legionella is known for its ability to form biofilms [34], and Brevundimonas has been associated with lead removal in wastewater treatment plants [35]. Anaerolineaceae are known for cadmium removal [36]. Several SRBs such as Desulfosporinus, Desulfovibrio, and Desulfobulbus, are known to proliferate in the presence of metals because they produce hydrogen sulfide that may precipitate with metallic ions [30]. Moreover, Clostridium seems to be a key player in this process by starting sulfate reduction and raising the pH of the environment in a context of metal contamination [37]. Most of the genera that were selected by metals (B TRGs), including those providing public goods, display species that were observed in wastewater treatment plants (WWTP) or in microbial communities in manure.
The resilience of the community would then have occurred through the early selection for metalresistant bacteria, providing common protection for the rest of the microbial community such as EPS formation or metal detoxification allowing cooperation and reducing the metallic stress; rather than the use of broad-host-range IncP plasmids. This is supported by the Q-PCR results. While the amount of the broadly distributed czcA gene was poorly affected by the experimental conditions, the pbrA gene, present within a narrower range of bacteria, was found to be enriched in metal conditions. In the TRG B1 and B4 groups, many Flavobacteriaceae and Brevundimonas species were observed and were previously found to influence Pb motility [35,49].
Microbial interaction: Linking bacteria acting in the public good with the rest of the community To assess the influence of bacteria acting as facilitator for the resilience of sedimentary microbiomes, we employed eLSA statistics pinpoint Hyphomicrobiaceae (Rhizobiales)as a keystone OTU that preceded or co-occurred with many other members of the microbiome. In the metal-contaminated microcosms, this OTU was replaced by another Rhizobiales, Xanthobacter of which several strains areknown to have an increased capacity to deal with metals50]. The central position of a Hyphomicrobiaceae and Xanthobacter OTUsin the microbiome may be explained by their capacity to fix nitrogen, which may beneficiate the rest of the community [51]. The subdivision 3 of Verrucomicrobia, previously identified in fresh-water sediments [52], was also found to be a key player of both control and experimental communities. In metal-contaminated sediments, Xanthobacter, Verrucomicrobia subdivision 3 member and other unidentified key OTUs displayed many links with public-good providing bacteria such as Zoogloea, Dechloromonas or SRBs. The latter facilitators may then be crucial for reducing the environmental pressure operated by metals for the keystone species of the community. In addition, Anaerolineaceae (TRG B3, Chloroflexi),, which was previously shown to be a key component of metal contaminated communities and was suggested to impact metal solubility [36], displayed a positive local similarity index with all phyla represented in the network revealing its importance for the rest of the microbial community. Anaerolineaceae cooccurred with, Methanobacterium (OTU 266) and Methanolinea (OTU 533), which are hydrogenotrophic methanogens Euryarchaeota. This interaction may then be highly beneficial for this group that was negatively impacted by metals from the beginning of the contamination. Finally, the interactions between OTU169 (TRG B4), member of Veillonellaceae, with the rest of the community, highlighted the risk of dispersion of antibiotic resistance genes as suggested in-situ [7] that might be co-selected with metal-resistance genes [53]. Indeed, Veillonellaceae are marker of fecal contamination and known to carry lots of antibiotic resistance genes [54], such as identified in Clostridia [10], or Anaerolineaceae found in WWTP [39].

Conclusion
The present microcosm experiment highlighted the presence of potential facilitator bacteria such as

Total and bioavailable metal content in sediments
For total metal quantification, 200 mg of dried and sieved (<63 μm) sediments were digested with concentrated hydrofluoric acid (Prolabo, Normapur) at 120°C for 48 h followed by aqua regia mineralization (hydrochloric + nitric acids, 2:1 v:v) at 120°C for 24 h. The total metal concentrations (Al, Fe, Cd, Cu, Pb, and Zn) were determined using inductively coupled plasma atomic emission spectroscopy (ICP-AES; Varian, Vista Pro, axial view) according to the concentration range [55].
A proxy of the bioavailable metals (Al, Fe, Cd, Cu, Pb, and Zn) fraction in the sediments was determined by a soft acid attack: 1 g of dried sediment was mixed with 20 mL of HCl solution (1 M) and agitated for 24 h. The samples were filtered with a cellulose acetate membrane (0.45 μm). The metal concentrations were measured using ICP-AES [55] and are shown in Supplemental Table S1.
High-throughput 16S rRNA gene sequencing Sequencing libraries were prepared using a dual-PCR setup, targeting variable regions V3 and V4 of the 16S rRNA gene, approx. 460bp. In the first step primers Uni341F (5'-CCTAYGGGRBGCASCAG-3') and Uni806R (5'-GGACTACNNGGGTATCTAAT-3') originally published by Yu et al. [56] and modified as described in Sundberg et al. [57] were used. In a second PCR reaction step the primers additionally included Illumina specific sequencing adapters and a unique combination of index for each sample [7,22]. the QIIME pipeline (https://github.com/maasha/qiime pipe) as previously described [22]. The MiSeq Controller Software was used to perform the sequence demultiplexing and sequencing adapters and primers were trimmed using cutadapt v1.18 [58] and any pairs were any of two primers couldn't be found were discarded. Trimmed reads were then processed using a custom BioDSL pipeline (http://maasha.github.io/BioDSL/). Specifically, paired-end reads were merged using merge_pair_seq, assembled pairs with a total length shorter than 300 bp or an average Phred score quality lower than Q25 were discarded. Clean reads were clustered in OTU using a 97% sequence similarity threshold using cluster_otus and USEARCH v7.0.1090 [59] and chimeric OTUs removed using UCHIME [60]. Each OTU cluster was annotated using mothur classify.seq() [61], an implementation of the RDP Classifier (method = wang) [62] against the Ribosomal Database Project Trainset 16 database [63].
Representative OTU sequence were aligned against a reference alignment using mothur and an approximate maximum likelihood phylogenetic tree was built using FastTree [64]. Samples (4 biological replicates x 2 technical replicates) with less than 3 000 OTU read counts were not included, as they do not provide enough coverage for further diversity analysis [65]. Information regarding the sequence counts for each sample is shown in Supporting Table S2.
Measuring the diversity through time Alpha-diversity estimates based on OTUs richness and the Shannon's diversity index were assessed using a rarefied dataset (R-package "Vegan" in RGui software [66]) containing 5150 sequences using the Past3 software [67]. The beta-diversity (a measure of species turnover over time) was estimated.
For this, a PERMANOVA test was performed using the R-package "Vegan" with a log 10 (1+x)transformed database (as OTU abundance variation was higher than 1 000-fold) with 10 000 permutations (Table S3). Differences in phyla between the control and metal-treated microcosms were assessed according to an FDR-corrected p-value inferred from a t-test (R-package "Vegan"),, and the different classes of Proteobacteria as well as the new order of β-Protobacteria were considered separately (Table S4). The evolution of phylogenetic relatedness in the control and metaltreated microcosms was then assessed between the in-situ profile and at each time point for both conditions using the phylogenetic tree obtained from OTU cluster representative sequence reads with the R-package "ape".. MPD quantifies the mean relatedness in a group of OTUs by considering all possible pairs of OTUs. MNTD quantifies the mean relatedness by considering only the nearest phylogenetic neighbors We calculated the weighted net relatedness index (NRI) and the nearest taxon index (NTI), which assess the number of standard deviations when comparing the mean pairwise distance (MPD) or the mean nearest taxon distance (MNTD), respectively, to the corresponding null distribution with 1 000 randomizations, 95% confidence interval, p < 0.05 [68,69] for each time point. An NRI or NTI value >2 indicates that the coexisting taxa are significantly closely related, while an NRI or NTI <-2 indicates significant phylogenetic overdispersion [70]. The metrics were calculated with the R-package "Picante" [71]. The between-community version of NRI (βNRI) and NTI (βNTI), βNTI [70], were determined (weighted, with 1 000 randomizations, R-package "MicEco" [72]) comparing each time point to time 0. A βNTI value >2 indicates a greater dispersion than expected by chance, while a βNTI <-2 indicates a static community. Intermediate values indicate that changes appear to be due to stochastic events [70].

Identification and validation of the time response groups
The 16S rRNA gene amplicon dataset was divided by time point, and the OTU profiles obtained after 0.5, 2.5, 3.5 and 6.5 months were assessed separately. The OTUs that responded significantly to metal stress were assessed using a negative binomial distribution and a generalized linear model, which were corrected with 1 000 resampling iterations of the residual variance (nbGLM, p < 0.05) [7].
The 600 selected OTUs were plotted in a generalized heatmap, and the Treatment Response Groups (TRG) were defined with a hierarchical cluster dendrogram (Euclidian distance and average clustering, R-package "gplots" [73]). The statistical validity of the obtained TRGs was tested against the null model by a Monte Carlo simulation including all OTUs [7,22].

Dynamic interaction network
The 16S rRNA gene amplicon dataset was then used to infer interaction networks between OTUs using an Extended Local Similarity Analysis (eLSA) [74]. The eLSA approach was applied using 1 technical replicate by biological replicate, to the most abundant OTUs in the 16S rRNA library with the eLSA implementation as previously described [75], using the following parameters: "lsa_compute data.txt result.txt -s 8 -r4 -p perm -x 1000 -d 4" [76][77][78]. The minimum LSA index was set at 0.8, and the pvalue threshold was set at 0.05. The network visualization and analysis was performed with Cytoscape 3.6.1 [79] using an organic-type layout and analyzed with the NetworkAnalyzer tool.
Functional and plasmid-associated gene content assessment by quantitative PCR.
The abundance of the metal resistance genes czcA and pbrA and the origin of transfer (oriT) of broad host range IncP plasmids were assessed by quantitative PCR. The primers used for the PCR amplification of czcA and oriT were retrieved from the literature (Roosa et al.,, 2014, Götz et al.,, 1996 (Table S5). The primers for the PCR amplification of pbrA were designed for this study. For pbrA, the genetic sequences of Cupriavidus metallidurans CH34, Klebsiella pneumonia and Enterobacter cloacae were aligned using the ClustalX. The best consensus sequences were used as templates for the design of degenerated primers. The physicochemical properties of each pair of primers were compared using the ThermoFisher Scientific Multiple Primer Analyzer tool [80]. The best candidates were selected (Table S5). The specificity of the primers was tested using BLAST (NCBI) and classical PCR of pure Cupriavidus metallidurans CH34 and Pseudomonas putida KT2440 cultures.