Located in the south of France, Toulon Bay is divided into two basins by an artificial seawall: the Little Bay, which is semi-enclosed, and the Large bay, which opens onto to the Mediterranean Sea (Supplementary Figure S1). The Little Bay is characterized by a high level of anthropic activities (Navy harbor, ferry transport, industries, wastewater sewage, and aquaculture). Due to its location, the Little Bay is directly impacted by pollution and is less affected than the Large Bay by offshore hydrodynamics. This confers to the Little bay a higher trace metal contamination than the Large Bay as demonstrated previously in the sediments and with punctual measurements in the water column (Tessier et al., 2011; Coclet et al., 2018).

Sampling was performed weekly from June 1st to 29th June, 2015. Two liters of seawater from the surface (1 m depth) and the bottom (1 m above sediment) of the water column were sampled at five sites across Toulon Bay (North-Western Mediterranean Sea, France), from the entrance of the bay to the north-western anthropized area (Supplementary Figure S1, Supplementary Table S1, and details in Supplementary Material). After the T1 sampling in site 41p, the French Navy performed an exercise that could have modified the water quality. In order to check that, we sampled again the water (T1b). The five sites were chosen to cover the whole range of dissolved metal contamination on the basis of trace metal distribution previously evaluated in Toulon Bay surface seawater (Coclet et al., 2018; Figure 1). Water samples were stored in a cooler and filtered through 0.2 μm polycarbonate membranes (Millipore). Filters were stored at −80°C until DNA extraction. Physico-chemical parameters including metals were weekly analyzed during the sampling period in order to characterize the environment (details in Supplementary Material). Standard oceanographic properties, including water temperature (°C), salinity, O₂ (mg L⁻¹ and %), pH, and chlorophyll a (chla) (μg L⁻¹) were measured on sites using a multi-parameter probe (Hydrolab DS5, OTT).

Subsamples of 10 mL of seawater were sampled at each site, depth, and date, filtered through a 90 μm nylon mesh, fixed with 0.25% (final concentration) glutaraldehyde on field and stored at −80°C until flow cytometry analysis. Autotrophic prokaryotes (Synechococcus-like), eukaryotic phytoplankton (picoeukaryotes and nanoeukaryotes) populations were characterized and enumerated using a BD Accuri^TM C6 (BD Biosciences) flow cytometer, as previously described (Coclet et al., 2018). Heterotrophic prokaryotes were enumerated after staining with SYBR Green as previously described (Cabrol et al., 2017).

DNA was extracted from the Millipore filters by a combination of enzymatic cell lysis (Ghiglione et al., 2009) and AllPrep DNA/RNA Mini Kit (QIAGEN) according to the manufacturer’s instructions. The protocol for the DNA extraction and the library preparation is fully described in Supplementary Material. We assessed bacterial community composition (BCC) by targeting the V4–V5 region of the 16S rRNA gene (Parada et al., 2016) and using Illumina Miseq 2 × 250 pb paired-end sequencing (Genoscreen, France).

Sequences were then demultiplexed and assigned to corresponding samples using CASAVA (Illumina). Forward and reverse reads were merged using PEAR 0.9.8 with default options (Zhang et al., 2014). Raw sequences were analyzed using MacQiime v.1.9.1 software (Caporaso et al., 2010). Briefly, barcode, primer, shorter sequences (<100 bp in length), and sequences with ambiguous base calls or homopolymer runs exceeding 10 bp were removed. The remaining sequences were assigned to operational taxonomic units (OTUs) and clustered at a 97% threshold using Uclust algorithm (Edgar, 2010), both closed and open reference OTU picking, based on the SILVA (release 128) database (Pruesse et al., 2007; Quast et al., 2013). Low abundance OTU (<0.005%) were filtered as recommended by Bokulich et al. (2013). Sequences classified as mitochondria or chloroplast were removed from the OTU table, corresponding to 3393 OTUs. A total of 349,009 reads were finally obtained representing 31,132 OTUs. OTU table was normalized by random subsampling to the smallest number of sequences (i.e., 4044). Samples with lower number of sequences were discarded (n = 16), 36 samples were remaining for following analyses (Supplementary Table S1). The 16S rRNA gene sequences have been deposited in the NCBI Sequence Read Archive (SRA) database under BioProject ID PRJNA514222¹.

All plots and statistical analysis were done in R RStudio (R Core Team, 2015). Alpha diversity calculations and error estimates were averaged, using QIIME script core_diversity_analyses.py (Caporaso et al., 2010), including estimates of rarefaction, Chao1, equitability, Shannon, and Simpson’s diversity. In order to compare the sites based on their chemical composition, taxonomic and functional community profiles, seawater samples were sorted in a non-metric multidimensional scaling (NMDS). Pairwise dissimilarities were calculated using Bray–Curtis metrics. Each resulting dissimilarity matrix was used to visualize sample differences via NMDS ordination using the “vegan” package in R (Oksanen et al., 2016). To test the null hypothesis that there were no significant difference between the groups discriminated according to sampling stations, sampling dates, and depth, similarities were analyzed by global and pairwise PERMANOVA tests, using “vegan” and “RVAideMemoire” packages, respectively.

To identify taxa that discriminated the sampling sites, the linear discriminant analysis (LDA) effect size method (LEfSe) (Segata et al., 2011) was used. This was performed with the LEfSe online tool in the Galaxy framework, using all default setting for data formatting and LDA effect size (Goecks et al., 2010). A similarity percentage (SIMPER) was run with a 90% cutoff and used to rank the percent contribution of individual biomarker to the dissimilarity between site differences (Clarke, 1993).

To investigate the relationships between BCC and measured environmental variables, Spearman rank correlation analysis, and redundancy analysis (RDA) were performed using the “vegan” package in R (Oksanen et al., 2016) (details in Supplementary Material). Briefly, for variable reduction and in order to create an efficient model from the most significant explanatory variables, vegan’s ordistep function were applied and among the 31 explanatory variables, 20 significant variables (P < 0.05) were kept for the following analyses The RDA model was tested by performing partial RDA, and variation partitioning to test the significance of the contribution of both groups of variables and each individual variable.

Interactions among the OTUs and between these OTUs and the 31 environmental variables were evaluated using CoNeT (Co-occurrence Network) plugin 1.1.1.beta (Faust et al., 2012), a plugin in a Cytoscape software 3.4.0. A similarity matrix was built with different metrics (Spearman correlation and Kullback–Leibler distance and a mutual information score) from OTUs. This initial network was re-defined by randomization. A permutation matrix representing a null distribution was obtained by resampling OTUs as described in Faust et al. (2012). In a permutation step, edge specific p-values were computed; however, for the final network, p-values of an edge were merged into one p-value according to Brown’s method. In the final step, the Benjamini–Hochberg multiple testing correction was applied (P < 0.05). Network characterization was evaluated using betweenness centrality (BC) and closeness centrality (CC) (details in Supplementary Material). Highly connected clusters were identified using the MCODE plugin version 1.4.0.beta2 (Smoot et al., 2011). Network characterization were evaluated using different topological indexes generated by Network Analyzer plugin.

Functional profiles were predicted from obtained 16S rRNA gene data using Tax4Fun (Aßhauer et al., 2015) based on KEGG category.

An increase in the trace metal concentrations was rationally observed from a reference site open on the sea with low metal concentrations (41p) to enclosed sites in the most anthropized area (MIS and 6ext) with high level of contaminations (Figure 1). The highest trace metal concentrations found in surface seawater were Zn, Cu, Pb, and Cd, while other metals/metalloids (As, Ba, Be, Cs, Cr, Mn, Sb, Sn, Ti, U, and V) were kept to concentration levels close to the geochemical background (Figure 1 and Supplementary Table S2). At the surface of seawater, Cd, Cu, Zn, and Pb concentrations were 6, 33, 35, and 48 times higher in 6ext water than in 41p water, respectively. These enrichment factors were quite similar (8, 25, 43, and 73 times) at the bottom of the seawater (Supplementary Figure S2). For both surface and bottom seawater, the concentrations of trace metals/metalloids were relatively stable during the sampling period (Figure 1, Supplementary Figure S2, and Supplementary Table S2). The most contaminated sites MIS and 6ext were significantly discriminated from the reference site 41p (PERMANOVA, P < 0.01) (Supplementary Figure S3A) and from sampling date (Supplementary Figure S3B) but not with depth (Supplementary Figure S3C).

Seawater temperature increased gradually during the studied period at all five sites, from 18.8°C–21.2°C in T0 to 22.0°C–24.3°C in TF (Supplementary Figure S4). Additionally, temperature was lower (1.9 ± 0.04°C) in bottom than in surface seawater. Surface salinity levels demonstrated low variations (38.4 ± 0.11). A very punctual drop was recorded in all sites at T2, down to 36.7 in MIS (Supplementary Figure S4). Chlorophyll a in surface seawater were lower in 41p (0.16 to 0.31 mg L⁻¹) than in other sites (0.59 to 2.08 mg L⁻¹). Additionally, chlorophyll a concentrations peaked in all sites at the middle of the survey. Except for 41p, chlorophyll a concentrations were higher at the bottom seawater than at the surface (Supplementary Figure S4). On average, in surface seawater, higher concentrations of DOC were found in 6ext (1.4 ± 0.07 mg L⁻¹) compared to 41p (1.2 ± 0.06 mg L⁻¹) (Supplementary Figure S4). Similarly, higher concentrations of TN were found in MIS (0.15 ± 0.02 mg L⁻¹) compared to 41p. (0.08 ± 0.02 mg L⁻¹) (Supplementary Figure S4). No significant difference was found in dissolved oxygen concentrations (mg L⁻¹, %), and pH between both all five sites and all sampling date (Supplementary Table S2). NO^3- and PO₄^3- concentrations of were below the detection limit (0.02 and 0.020 μM, respectively) in a majority of samples. As a consequence, they are not presented and will not be considered further.

Heterotrophic bacterioplankton abundances estimated using flow cytometry are shown in Supplementary Figure S5. We observed an increase in the bacterioplankton abundance from the reference site 41p (3.8 × 10⁵ to 4.2 × 10⁵ cell mL⁻¹) to the most enclosed site 6ext (8.0 × 10⁵ to 1.4 × 10⁶ cell mL⁻¹) in both surface and bottom seawater. Additionally, heterotrophic bacterioplankton abundance was stable in 41p during the experiment whatever the depth, while strong dynamic was observed in both surface and bottom 6ext seawater, with highest abundances in T2 and T3, respectively.

Over the sampling, observed richness, and Chao1 index (Supplementary Figure S6) were significantly higher in the uncontaminated site 41p and declined to less than half compared to 6ext (ANOVA, P < 0.01). In contrast, Shannon, Simpson, and evenness indexes did not significantly vary between sampling sites (Supplementary Figure S6).

NMDS results showed a significant segregation of bacterioplankton communities between both sampling sites (PERMANOVA; P = 0.001) and dates (PERMANOVA; P = 0.024) but not with depths (PERMANOVA; P = 0.983) (Figure 2A and Supplementary Table S3). More precisely, pairwise PERMANOVA tests (Supplementary Table S3) confirmed significant differences between the reference 41p site and all other sites (PERMANOVA; P < 0.05) and significant differences between 6ext and Pt12–Pt15 (PERMANOVA; P < 0.05).

The most abundant phyla ( > 1% of all sequences across all samples) were Proteobacteria (53 ± 1.1%), followed by Bacteroidetes (28 ± 1.1%), Cyanobacteria (9.1 ± 0.7), Actinobacteria (6.1 ± 0.6), and Verrucomicrobia (3.0 ± 0.3%). At the family level (Figure 2B), Flavobacteriacea comprised the majority (22 ± 1.1%) of the total sequences, followed by Alphaproteobacteria groups SAR11 Surface 1 clade (18 ± 1.1%), Rhodobacteraceae (13 ± 0.65%) and SAR116 clade (6.1 ± 0.49%). Cyanobacteria, Actinobacteria, and Verrucomicrobia phyla were dominated by Synechococcus Family I (9.1 ± 0.72%), Microbacteriaceae (4.9 ± 0.71%), and Puniceicoccaceae (1.7 ± 0.19%). Archaea represented less than 0.1% of the total sequences, thus they were considered as “Others” in this study and variations of this group have not been taken into consideration.

Clear shifts in the relative abundance of individual bacterioplankton OTUs were observed between sampling sites. Using LEfSe analysis, we determined 43 differentially abundant taxa between four sites. 41p and 6ext, the two most geochemically contrasted sites, exhibited the highest number of taxonomic indicators (25 and 6, respectively) (Figure 3). For 41p site, the most significantly enriched bacterioplankton sequences were Acidimicrobiia Candidatus Actinomarina, members of the family SAR406, 7 Alphaproteobacteria OTUs and Gammaproteobacteria Thiothrix. Actinobacteria Candidatus Aquilina and both Flavobacteriia Formosa and Betaproteobacteria Hydrogenophilaceae were found to be overrepresented in MIS and 6ext, respectively. SIMPER analysis revealed that, on average, 13 biomarkers were identified as major contributors to differences between sites detected in the NMDS. The higher number of biomarkers identified as major contributors was found between 41p and the other sampling sites (between 17 and 32 biomarkers) (Supplementary Table S4). Overall, each biomarker contributed to at least 1% of the dissimilarity between sites. Candidatus Aquilina, PS1 clade, NS7 marine group, Planktomarina and Tenacibaculum, notably, drove together more to 32% of disparity between cluster 41p and the other clusters (Pt15, 6ext, and MIS).

Among the 31 environmental variables available for multivariate analysis, 20 parameters were selected by pairwise correlations. The first and the second axes of RDA captured 47.1% and 6.4% of BCC variance, respectively (Figure 4A).

Variation partitioning provided more details about the relative contributions of nutrients, contaminants, terrigenous, and “others” factors in the observed spatial changes in BCC (Figure 4B). The combination of the four sets of explanatory variables explained 41% of BCC variation across sites. Among the set of explanatory variables, contaminants were the main contributors to variations in BCC (30%). Contaminants also appeared to have significant influence on BCC through their interaction with the other sets of variables. The contribution of terrigenous (19.6%), nutrient and biotic (15.0%) and “marine tracers” (21.6%) variables were also significant but to a lower extent than the contribution of contaminants variables. The highest individual contributions were attributed to Mn (1.2%), DOC (1.2%), salinity (1.1%), Zn (1.1%), and Cd (1.0%) (Supplementary Table S5). Thus, among the five explanatory variables which contributed the most to the biological variation, three variables were contaminants (Mn, Cu, and Cd).

Network analysis resulted in a global network representing 158 OTUs and 21 environmental variables significantly connected. Among 1664 edges, 92% were inter taxa edges, and 8% represented taxon-environment interactions. Most of the connections (67%) corresponded to negative interactions. To visualize the results, eight subnetworks were extracted from the global network (Figure 5 and Supplementary Figure S7).

Globally, by plotting BC vs. CC, 12 OTUs with the highest values for both parameters (arbitrarily determined as BC > 0.02 and CC > 0.5) were selected, representing putative keystone species within Toulon Bay (Supplementary Table S6). Both of these indices indicate that SAR11 Surface 1 clade and Rhodobacteraceae (Alphaproteobacteria), Synechococcus (Cyanobacteria), Candidatus Aquiluna (Actinobacteria), and Balneola (Bacteroidetes) were central in this network. Most of Alphaproteobacteria, including SAR11, SAR116, and Rhodobacteraceae, and Cyanobacteria Synechococcus were positively correlated with each other, but negatively correlated with Bacteroidetes (NS11-12 and NS5 marine group, Cryomorphaceae and Balneola), and Actinobacteria (Candidatus Aquilina).

In the subnetworks, there were 19 correlations with contaminants (8 with Pb, 6 with Cu, 3 with Cd, and 2 with Zn) and two correlations with chlorophyll a. Several taxa that were identified as being key drivers of community dissimilarity using LEfSe analysis were highly correlated to contaminants. In particular, the Alphaproteobacteria Rhodobacteraceae and SAR11, the Bacteroidetes Balneola and the Cyanobacteria Synechococcus showed negative correlations with at least one of the contaminants (Cd, Cu, Pb, or Zn). Conversely, Actinobacteria Candidatus aquilina and Bacteroidetes NS11-12 marine group were positively correlated with Pb and Cu, respectively. Gammaproteobacteria Haliaceae was positively correlated with Cd.

The functional profiles of bacterioplankton communities was predicted using Tax4Fun (Aßhauer et al., 2015). Predicted functions were classified as KEGG orthologs (KOs) resulting in the identification of 3619 KOs across all samples, using on average, 60 ± 12% of total OTUs (Supplementary Table S7). Clustering showed that the differentiation of the sampling sites based on microbial predicted function profiles followed the chemical gradient described earlier, from the reference site toward the most enclosed and contaminated 6ext site (Figure 6). Compared to the uncontaminated site (41p), used here as a reference, some KEGG pathways seemed overrepresented in contaminated sites: Carbohydrate metabolism, Xenobiotics biodegradation and metabolism, Metabolism of terpenoids, and polyketides. Other KEGG pathways were found to be underrepresented as “Cell growth and death,” “Nucleotide metabolism,” “Energy metabolism,” and “Metabolism of cofactors and vitamins.”

Marine coastal areas are seen as valuable sentinel ecosystems, providing signals that reflect anthropogenic chemical pollution (The MerMex Group et al., 2011). Concentrations of trace metals in the Toulon Bay were found to be at least comparable or higher to those reported in aquatic ecosystems receiving high contamination loads (Boyd et al., 2005; Pringault et al., 2012). The high enrichment factors (6- to 234-fold) compared to the geochemical background of the Mediterranean Sea (Morley et al., 1997) indicated significant anthropogenic inputs that could be attributed to numerous recent (e.g., large boat traffic, harbor activities, antifouling coatings) (Turner, 2010) and historical events (2nd World War) (Tessier et al., 2011). Thus, we hypothesized a potential effect of trace metals on the bacterioplankton communities.

As previously reported in coastal waters (Nayar et al., 2004; Caroppo et al., 2006; Pringault et al., 2016), bacterioplankton was more abundant in the most contaminated zones (MIS and 6ext) than in the uncontaminated site 41p. Bacteria have been reported to have a wide range of metal detoxification or resistance mechanisms by sequestering, excluding or precipitating metals (Nies, 1999). The higher abundance in contaminated sites could therefore suggest that the community is well adapted. Moreover, bioavailability of trace metals could also contribute to explain the absence of bacterioplankton abundance decrease in the most contaminated sites in term of abundance. Trace metals bioavailability is strongly controlled by the presence of dissolved organic ligand, especially in marine waters (Louis et al., 2009; Bruland et al., 2013). We considered here DOC concentrations and the dissolved fraction of trace metals, which is closer to the bioavailable fraction than the total concentration (Paquin et al., 2002; Cindrić et al., 2017) and consequently allow a more reliable assessment of their roles. Nevertheless, we neither evaluated the speciation of the main trace metals nor DOC composition in this study. More work will be thus needed to state about the respective contributions of bioavailability modulation and community adaptation in the bacterial abundance increase in Toulon Bay.

While bacterioplankton diversity did not show any significant trend along the contamination gradient, the community richness (Chao1 index) was shown to decrease progressively from the uncontaminated site (41p) toward the most contaminated sites (MIS and 6ext). A similar pattern was observed in South-Eastern Australia in metal- and PAH-contaminated sediments (Sun et al., 2012) as well as in contaminated sediments in the Northern Adriatic Sea (Korlević et al., 2015). This may be interpreted in light of classis disturbance theory that predict that in highly disturbed ecosystems, richness would decrease and a few species would dominate (Odum, 1985; Allison and Martiny, 2008). As a matter of fact, sensitive species are expected to be lost whereas tolerant species benefit from emptying ecological niches (Atlas et al., 2015). Thus, the multi-contamination appeared to have disruptive effect on the alpha diversity of the bacterioplankton communities.

Gillan et al. (2005) demonstrated that alpha-diversity alone is a poor indicator of ecosystem stress in chronically polluted systems, as the proliferation of new tolerant species can lead to a recovery of the diversity. As in many coastal sites, microbial assemblages in our study were dominated by a few families Flavobacteriaceae, SAR11 Surface 1 clade, Rhodobacteraceae, SAR116 clade, Synechococcus Family I, and Microbacteriaceae. This is a core of generalist families in marine seawater and its members have been found abundant in numerous other marine coastal studies (Chow et al., 2013; Cram et al., 2015). Despite overall similarities, BCC was strongly affect by trace metal contamination, the five sampling sites clearly exhibited distinct microbial community’s structures, along low-to-high contamination continuum. This observation is also in agreement with the disturbance hypothesis, as discussed for alpha diversity, and tends to complement previous observations of the structural role of contaminants for bacterioplankton in other coastal areas (Pringault et al., 2016; Marisol et al., 2018; Wang et al., 2018).

The uncontaminated site (41p) was characterized by higher relative abundances of OM1 clade (Acidimicrobiia), SAR11, Rickettsiales, Rhizobiales, and Rhodobacteraceae (Alphaproteobacteria), Thiothrichaceae (Gammaproteobacteria), and Synechococcus (Cyanobacteria). Correspondingly, a number of previous works reported that Synechococcus were sensitive to high trace metal levels (Cassier-Chauvat and Chauvat, 2014; Coclet et al., 2018). Additionally, SAR11, Rickettsiales, Rhizobiales, and Synechococcus are well known to be characteristic of marine zones away from anthropogenic influence (Gilbert et al., 2009, 2012; Fuhrman et al., 2015).

The most contaminated sites (MIS and 6ext) were enriched in Candidatus Aquiluna (Actinobacteria), Formosa and Tenacibaculum (Flavobacteriia) as well as one from Sphingobacteriales NS11-12 marine group (Bacteroidetes). Networks analysis showed also that Planktomarina (Alphaproteobacteria), Hydrogenophilaceae (Betaproteobacteria) and Halieaceae (Gammaproteobacteria) were positively associated with contaminants, in agreement with the distribution of these groups along the contamination gradient. Previous studies in freshwater or marine sediment reported that Hydrogenophilaceae and Candidatus Aquiluna (Microbacteriaceae) were positively correlated to high trace contaminant levels (Acosta-González et al., 2013; Garris et al., 2018). The increase of Actinobacteria and Betaproteobacteria in the most contaminated site was not unexpected, as they are widely recognized to play a key role in oil degradation (Acosta-González et al., 2013) and heavy metals transformation, especially, in the case of nutrient enrichments (Garris et al., 2018). Microbacteriaceae (Actinobacteria) were widely reported as members of highly copper-polluted sediment (Besaury et al., 2014) or Zn-polluted freshwater (Ni et al., 2016). Albarracín et al. (2010) have shown that Actinobacteria members isolated from copper-contaminated sediment are able to significantly diminish the bioavailable copper throughout bioaccumulation. Members of the Bacteroidetes members, especially Flavobacteria families are also known to represent one of the most abundant groups of bacteria in coastal areas, especially in contaminated sites due to their ability to tolerate toxic effects of certain metals (Sun et al., 2013: Korlević et al., 2015). The Bacteroidetes NS11-12 marine group have been detected mainly in marine habitats (Meziti et al., 2015), unfortunately without clear ecological implications. At a finer scale, for most of clusters, members of both Alphaproteobacteria (mainly SAR11) and Bacteroidetes (mainly Flavobacteriia) groups describe above, exhibit the highest number of edges with high level of connections within each group in the network analysis. These important findings display that Alphaproteobacteria would highly influence the functioning of low contaminated site while Bacteroidetes would highly influence the functioning of contaminated zones.

Most of studies on bacterioplankton communities have focused on community diversity and composition responses to perturbation (Jeffries et al., 2016). The effect of chemical contaminations on bacterial community functionalities have only been investigated in freshwater ecosystems (Wang et al., 2018; Sjöstedt et al., 2018). In a highly polluted area such as Toulon Bay, the functional stability of microbial communities could be challenged. Major predicted bacterial functions including those necessary to basic metabolisms were conserved across the chemical gradient. Functional recovery occurred because of species replacement or change in relative abundance of taxa, i.e., functionally redundant rare species could become abundant in response to the contaminations. This is in agreement with the replacement scenario (Comte and Del Giorgio, 2011), and previous findings showing that the rare biosphere is an important reservoir for recruitment of bacterial taxa during environmental change (Sjöstedt et al., 2012; Comte et al., 2014). However, both a decrease in some critical functions (cell growth and death or energy metabolism) and an increase in functions associated to multiple stress resistance (e.g., xenobiotic biodegradation and metabolism, signaling molecules and interactions, secondary metabolites metabolism) remarkably exhibited that bacterial communities in the most polluted sites developed specific functions to face their environment. As the sites are not physically separated, there was no discontinuity in both communities and functions but a clear gradient of the relative components throughout the metal gradient. Metal contaminated environments were shown to exert a high selective pressure toward the transfer of several genes, developing thus resistance against high metal concentrations (Azarbad et al., 2016). Only two field studies on long-term metal polluted areas have shown shifts in microbial community structure and functions or increased occurrence of resistance genes that have made microbial communities resistant to toxic metals concentrations (Hemme et al., 2010; Kang et al., 2013). However, both taxonomic affiliation of metal-specific metabolic traits and the currently available KEGG database are insufficient to validate the hypothesis that heavy metal stress is a key environmental factor shaping the function of microbial communities.

Although trace metal gradient seemed to be the main driver of the bacterioplankton community features, other local selective pressures could explain a significant proportion of their variance. Indeed, bacterioplankton communities in coastal areas greatly vary in space because of sharp gradients in salinity, nutrients, depth, among other properties (Fuhrman et al., 2006; Gilbert et al., 2009; Fortunato et al., 2012).

A very similar trend could have been observed if DOM quantity or quality varied and became more labile in the most contaminated sites. Indeed, an increase in DOM lability often leads to the development of a copiotrophic response of the bacterial community, selecting a limited number of bacterial taxa, and temporarily reducing microbial diversity (Nelson and Wear, 2014; Pedler et al., 2014). The variation partitioning performed in our study is in agreement with this hypothesis, since DOC appeared as one of the most influent drivers of community structure variations. Thus, organic matter spatio-temporal dynamics could indeed be an important driver of microbial communities in Toulon Bay. However, according to the literature, only a part of the bacterial taxa significantly enriched in the most contaminated sites of Toulon Bay is known as copiotroph. Indeed, while OM60 (NOR5) clade (Halieaceae), Candidatus Aquiluna (Microbacteriaceae), and Planktomarina (Rhodobacterales) are known to benefit from labile DOM (Yan et al., 2009; Spring et al., 2015), Alphaproteobacteria members (SAR11 and AEGEAN-169 marine group) and Bacteroidetes (Balneola and Aquibacter) are rather known as oligotrophic groups (Cram et al., 2015; Sipler et al., 2017). Copiotrophs are expected to outcompete oligotrophs, so the fact that presence of groups identified as copiotrophs suggest that each has capitalized on reduced competition following trace metal contamination. Taken together, all these arguments point out the major and predominant influence of the multi-contamination gradient on microbial diversity and represent a significant demonstration of the potential influence of human activities on marine microbial life. Since we only considered inorganic contaminants, we cannot rule out the potential influence of organic contaminants such as PAH, PCBS, or organometals which were already measured in the sediment of Toulon Bay (Misson et al., 2016).

Organic contaminants as polyaromatic hydrocarbons could explained more precisely the distribution of OTUs (Qian et al., 2017). This hypothesis could be confirmed by the strong enrichment in xenobiotics biodegradation and metabolism, carbohydrate metabolism and metabolism of terpenoids and polyketides in the highest contaminated sites. Naphthalene degradation, benzoate, and aminobenzoate biodegradation were found to be higher in the most contaminated zones and is consistent with the finding of previous studies (Wang et al., 2018).

Finally, another striking result is that the multi-contamination gradient in the water column of Toulon Bay represented a significant perturbation, affecting microbial diversity all along our survey. This observation differs from what was observed previously in surface sediment of Toulon Bay (Misson et al., 2016), suggesting that the contemporary contamination in the water affects more strongly microbial communities that the large historical contamination trapped in the sediment. Such difference could be linked to the probably higher dispersal rate in the water column than in the sediment, i.e., 3.4 days for a total renewal of the water of Toulon Bay (Dufresne et al., 2014). Thus, this higher dispersal rate could tend to reduce the adaptive abilities of bacterioplankton, leading to a higher sensitivity to contamination when compared to the longer-term stability of the sediment.

The network indices (BC and CC) and large numbers of nodes (179) and edges (1664) between taxa in our network analysis suggest that SAR11, Rhodobacteraceae, Synechococcus, and Microbacteriaceae play a key role in Toulon Bay seawater. Despite some connections between OTUs and trace metals, there were far more numerous connections among OTUs. This suggests that trace metals significantly influence the dynamics of only a few microbial groups, and could rather influence indirectly, via biotic interactions, the whole community in Toulon Bay seawater. Community structure is known to be influenced by interaction among species (Lima-Mendez et al., 2015). Some of these interactions could be cooperative interspecific interactions (Eiler et al., 2012; Chow et al., 2013) or competitive interactions (i.e., competition, niche partitioning, grazing, or parasitism) (Fuhrman et al., 2008; Berdjeb et al., 2018) in a multispecies community.

Keystone OTUs play important roles in the regulation of network interactions (Wright et al., 2012) and their loss may increase a community’s vulnerability to perturbation. Thus, the strong decline of Alphaproteobacterial and Cyanobacterial keystone OTUs, especially, along the trace metal contamination gradient could have indirect effect on the BCC by modifying biotic interactions and favor opportunist or tolerant taxa (Lawes et al., 2016). This is in agreement with community turnover and the development of copiotrophic groups highlighted above.

Additionally, in network analysis, positive correlation between on the one hand pico- and nanophytoplankton abundances, as well as chlorophyll a, and on the other hand trace metals concentrations describe potential conditions that may favor specific groups of bacterioplankton. Indeed, Fuhrman et al. (2015) explained that weekly timescale is appropriate for studying bacterial dynamics associated with phytoplankton, which have a direct influence on the composition of bacterial communities via cross-feeding interactions or the effect of toxins, as well as indirect influence (e.g., oxygen depletion). Phytoplankton succession may influence the dynamics of the bacterioplankton community, as differences in the phytoplankton structure lead to differences in quantity and quality of exuded organic matter that can be used by the bacterial community for growth (Marisol et al., 2018). Consequently, effects of pollutants on phytoplankton or bacterioplankton might have indirect consequences for the counterpart, depending on the possible interactions between both compartments.

In summary, Toulon Bay offered a remarkable study area allowing to address trace metal multi-contamination impacts on prokaryotic community shaping throughout contrasted but connected sites. Where trace metal contamination increased, bacterioplankton abundance was higher together with a reduction of only community richness. Seawater contamination induced significant bacterioplankton communities shifts, with Alphaproteobacteria (SAR11, Rickettsiales, Rhizobiales, and Rhodobacteraceae) and Synechococcus dominance where trace metal concentrations were the lowest, whereas Flavobacteria, Betaproteobacteria (Hydrogenophilaceae) and Gammaproteobacteria (Halieaceae), especially including copiotrophic organisms, dominated in highly contaminated sites. In addition, this study showed that human-induced disturbances can have an effect on predicted functional potential of marine coastal bacterioplankton, and that the response can differ according to the type of functions. However, observed functional redundancy lead to think that the effect of trace metals could be complex, highly context-dependent and depending on biotic interactions. While our results provide a better understanding of anthropogenic influence on coastal ecosystems and the response of microbial communities, future studies are needed to elucidate the contribution of DOC, phytoplankton-bacteria specific interactions and organic pollutants on bacterioplankton communities in Toulon Bay.