Sampling methods and microcosm maintenance
Sediments were sampled in the Senséecanal (Férin, France, 50°18’39.0’‘N—03°05’05.4’‘E) on the fifth of December 2016 using a plexiglass core sampler (Ø = 7 cm) to avoid the mixing of sediment layers. Thirty sediment cores were collected, and the top first 1.5 cm were retrieved, homogenized and stored on ice during transport. Freshwater was sampled and stored at 4°C. On the sampling day, eight polypropylene microcosms (750 cm³) were filled with 120 mL of sediment and 500 mL of freshwater. Every day, 450 cm³ of the water column of each microcosm were replaced by water sampled in the river monthly, stored and oxygenized at 4°C with no sterilization step. During the first eight weeks of maintenance, a mixture of metals (ZnCl2, CuCl2(2H2O), PbCl2 and CdCl2(1/2H2O) was gradually added to four metal-treated microcosms at a rate of 20 mL every four days. Final metal concentrations are displayed in table S1. The pH remained stable over the course of the experiment (ca. pH 8.0). Microcosms were maintained in an air-conditioned room (23°C) over 6.5 months. The sediments were sampled for DNA extraction after two weeks (T0.5) and every month thereafter (T1.5, T2.5, T3.5, T4.5, T5.5 and T6.5).
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 .
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  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.  and modified as described in Sundberg et al.  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]. First PCR reaction used Phusion® High-Fidelity DNA Polymerase (Thermo scientific, Carlsbad, CA, USA) with an initial denaturation step at 98 °C for 30 s, followed by 35 amplification cycles consisting of denaturation at 98 °C for 10 s followed by 30 s at 61 °C and 15s at 72 °C and final extension of 7 min at 72 °C. Second PCR reaction used PCRBIO HiFi Polymerase (PCR Biosystems Ltd, UK) with an initial denaturation step at 95 °C for 1 min, followed by 15 amplification cycles consisting of denaturation at 95 °C for 15 s followed by 15 s at 56 °C and 30 at 72 °. After each PCR reaction, amplicon products less than 200 bp were purified using HighPrep™ PCR Clean-up System magnetic beads (MagBio, MD, USA) according to the manufacturer’s instructions. The samples were adjusted to equimolar concentrations using SequalPrep Normalization Plate 96-well kit (Thermo Fisher Scientific). A final library pool was made using 5µL of each individual normalized sample library, concentrated using DNA Clean TM and Concentrator 5 kit (Zymo Research, Irvine, CA, USA) in a final volume of 10µL. Finally, the pooled libraries were adjusted to 4nM and sequenced using MiSeq v2 reagents kit (500 cycles) using 2 × 250 bp paired-end mode on an Illumina MiSeq platform (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. The amplicon sequences were analyzed with the QIIME pipeline (https://github.com/maasha/qiime pipe) as previously described . The MiSeq Controller Software was used to perform the sequence demultiplexing and sequencing adapters and primers were trimmed using cutadapt v1.18  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  and chimeric OTUs removed using UCHIME . Each OTU cluster was annotated using mothur classify.seq() , an implementation of the RDP Classifier (method = wang)  against the Ribosomal Database Project Trainset 16 database . Representative OTU sequence were aligned against a reference alignment using mothur and an approximate maximum likelihood phylogenetic tree was built using FastTree . 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 . 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 ) containing 5150 sequences using the Past3 software . 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 log10(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 metal-treated 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 . The metrics were calculated with the R-package "Picante” . The between-community version of NRI (βNRI) and NTI (βNTI), βNTI , were determined (weighted, with 1 000 randomizations, R-package "MicEco” ) 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 .
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) . 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” ). 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) . 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 , using the following parameters: “lsa_compute data.txt result.txt -s 8 -r4 -p perm -x 1000 -d 4” [76–78]. The minimum LSA index was set at 0.8, and the p-value threshold was set at 0.05. The network visualization and analysis was performed with Cytoscape 3.6.1  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 . 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.
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, czcA and pbrA. 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, czcA or pbrA) 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) 
Ratio = (Etarget)ΔCP target (T0—Tx) / (Eref)ΔCP ref (T0—Tx) (1)
where E is the PCR efficiency calculated according to Ramakers et al.,, 2003, and CP is the point at which the amplification curve crosses the threshold.
Duplicate PCRs were performed for each DNA sample and target gene combination using SYBRTM Green master mix (Applied BiosystemsTM). Each replicate consisted of serial dilutions (0.5, 0.25, 0.125 and 0.0625 ng/µL for oriT and pbrA and 1, 0.7, 0.5, 0.25 ng/µL for czcA). After an initial denaturation step at 94°C for 10 min, 40 cycles of 15 seconds at 95°C, 30 seconds at 60°C and 30 seconds at 72°C were performed followed by a final denaturation step (60 to 90°C, +0.3°C/min).