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).

# Bacterial Community Composition (BCC) in Relation With Environmental Variables

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.

FIGURE 2 | Spatial variations of bacterioplankton community structure. (A) Non-metric dimensional scaling (nMDS) ordination based on Bray–Curtis dissimilarity for 16S rRNA gene libraries between the different sampling sites. Each symbol corresponds to a distinct site, time, and depth point. (B) Bar plot showing the average family-level contribution for each site. Families that on average comprised less than 1% of the libraries are grouped as "Others." Colors indicate different classes within phyla.

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.

### Bacterioplankton Predicted Functional Capabilities of Microbial Communities

The functional profiles of bacterioplankton communities was predicted using Tax4Fun (Aßhauer et al., 2015).

FIGURE 4 | Contributions of environmental variables to spatial differentiation in bacterioplankton community structure. (A) Redundancy analysis (RDA) ordination diagram of the first two axes for bacterioplankton community composition. The percentage of the spatial variation in community structure explained by each axis is indicated in parentheses after the axis label. The constrained sets of environmental variables analyzed are indicated as vectors. NanoEuk: Nanoeukaryote abundance; PicoEuk: Picoeukaryotes abundance. (B) Venn diagram showing the variation partitioning (%) of the spatial variations in bacterioplankton community structure among four environmental datasets (Contaminants, Terrigenous elements, Marine tracers, and Nutrient and Biotic).

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."

 This study was performed in a multi contaminated coastal area, Toulon Bay, which presents far wider contamination gradients than natural environmental gradients. By combining seawater physicochemical and geochemical analyses with patterns of bacterioplankton diversity and predicted functions, we provide an innovative insight into bacterioplankton ecology allowing to better understand how trace metals could shape communities.

# Trace Metal Contaminations as Drivers of Bacterioplankton Communities

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; Cindric 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 metaland PAH-contaminated sediments (Sun et al., 2012) as well as in contaminated sediments in the Northern Adriatic Sea (Korlevic 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 lowto-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: Korlevic 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.

### Implication Regarding Others Local Selective Pressures

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.

# Biotic Interactions as Evidences of Indirect Effects of Trace Metals

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.

### ORIGINALITY-SIGNIFICANCE STATEMENT

Human-induced chemical contamination gradients challenge our understanding of microbial communities' ecology, especially in the very dynamic marine coastal environment. Through the spatio-temporal investigation of microbial diversity dynamics in a model ecosystem, this work highlights contributions of anthropogenic trace metal contamination to marine bacterioplankton ecology. Identify the key aspects of originality and significance that place the work within the top 10% of current research in environmental microbiology.

### AUTHOR CONTRIBUTIONS

CG, J-FB, and BM proposed and designed the study. CC, CG, GD, CLP, J-UM, J-FB, and BM organized and performed field sampling. CC, CG, GD, DO, SD, CLP, and BM analyzed the samples. CC, CG, J-FB, and BM interpreted the results. All authors participated in manuscript redaction, thus approving the publication of the content.

### FUNDING

This study was supported by PREVENT research program [funded by Toulon University, Toulon-Provence-Méditerranée (TPM) and the Conseil Départemental du Var]. CC was funded by the Provence-Alpes-Côte d'Azur region and supported by TPM.

# ACKNOWLEDGMENTS

The authors wish to thank the French Navy and the LASEM for diving and sampling assistance, Dr. G. Culioli (MAPIEM, Toulon, France) and L. Favre (MAPIEM, Toulon, France) for sampling and sample treatment assistance.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2019. 00257/full#supplementary-material

### REFERENCES

fmicb-10-00257 February 20, 2019 Time: 17:14 # 12

environmental gradients in the Red Sea. ISME J. 11, 138–151. doi: 10.1038/ ismej.2016.99

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Coclet, Garnier, Durrieu, Omanovi´c, D'Onofrio, Le Poupon, Mullot, Briand and Misson. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Shift in Natural Groundwater Bacterial Community Structure Due to Zero-Valent Iron Nanoparticles (nZVI)

Marc Crampon\*, Catherine Joulian, Patrick Ollivier, Mickaël Charron and Jennifer Hellal

Bureau de Recherches Géologiques et Minières, Orléans, France

Toxic and persistent contaminants in groundwater are technologically difficult to remediate. Remediation techniques using nanoparticles (NPs) such as nZVI (Zero-Valent Iron) are applicable as in situ reduction or oxidation agents and give promising results for groundwater treatment. However, these NP may also represent an additional contamination in groundwater. The aims of this study are to assess the impact of nZVI on the nitrate-reducing potential, the abundance and the structure of a planktonic nitrate-reducing bacterial community sampled in groundwater from a multicontaminated site. An active nitrate-reducing bacterial community was obtained from groundwater samples, and inoculated into batch reactors containing a carbon substrate, nitrate and a range of nZVI concentrations (from 0 to 70.1 mg Fe.L−<sup>1</sup> ). Physical (pH, redox potential), chemical (NO<sup>−</sup> 3 concentrations) and biological (DNA, RNA) parameters were monitored during 1 week, as well as nZVI size distribution and mortality of bacteria. Nitrate-reducing activity was temporally stopped in the presence of nZVI at concentrations higher than 30 mg L−<sup>1</sup> , and bacterial molecular parameters all decreased before resuming to initial values 48 h after nZVI addition. Bacterial community composition was also modified in all cultures exposed to nZVI as shown by CE-SSCP fingerprints. Surprisingly, it appeared overall that bacteria viability was lower for lower nZVI concentrations. This is possibly due to the presence of larger, less reactive NP aggregates for higher nZVI concentrations, which inhibit bacterial activity but could limit cell mortality. After 1 week, the bacterial cultures were transplanted into fresh media without nZVI, to assess their resilience in terms of activity. A lag-phase, corresponding to an adaptation phase of the community, was observed during this step before nitrate reduction reiterated, demonstrating the community's resilience. The induction by nZVI of modifications in the bacterial community composition and thus in its metabolic potentials, if also occurring on site, could affect groundwater functioning on the long term following nZVI application. Further work dedicated to the study of nZVI impact on bacterial community directly on site is needed to assess a potential impact on groundwater functioning following nZVI application.

Keywords: zero-valent iron nanoparticles (nZVI), groundwater, nitrate-reducing bacterial community, ecotoxicity, remediation

### Edited by:

Fabrice Martin-Laurent, Institut National de la Recherche Agronomique (INRA), France

### Reviewed by:

Naresh Singhal, The University of Auckland, New Zealand Ana Isabel Pelaez, University of Oviedo, Spain

> \*Correspondence: Marc Crampon<EMAIL_ADDRESS>
### Specialty section:

This article was submitted to Microbiotechnology, Ecotoxicology and Bioremediation, a section of the journal Frontiers in Microbiology

Received: 15 September 2018 Accepted: 01 March 2019 Published: 19 March 2019

### Citation:

Crampon M, Joulian C, Ollivier P, Charron M and Hellal J (2019) Shift in Natural Groundwater Bacterial Community Structure Due to Zero-Valent Iron Nanoparticles (nZVI). Front. Microbiol. 10:533. doi: 10.3389/fmicb.2019.00533

# INTRODUCTION

fmicb-10-00533 March 19, 2019 Time: 11:44 # 2

Due to their specific properties and high specific surface area, engineered nanoparticles (NP) are increasingly produced and used, leading to higher releases in the environment (Buzea et al., 2007). In particular, use of nanotechnology for environmental applications, such as remediation of contaminated groundwater, have steadily increased (Bardos et al., 2018). Among NP, nanoscale Zero-valent Iron (nZVI) particles are one of the most widely used NP for nanoremediation, because of their capacity to degrade a wide range of contaminants (Zhang, 2003; Zhang and Elliott, 2006; Mueller et al., 2012). Zerovalent iron (ZVI) is efficient for the dechlorination of PCE because of its ability to dehalogenate chlorinated compounds by chemical reduction (Fu et al., 2014). ZVI can be used in groundwater in permeable reactive barriers (PRB), that show effective degradation, but this technique is expensive, invasive and constrained by installation limits (Obiri-Nyarko et al., 2014). In this framework, an interesting alternative, can be direct injection of ZVI NPs into the aquifer, enabling the treatment of both source and plume areas (Cundy et al., 2008; Thiruvenkatachari et al., 2008; Stefaniuk et al., 2016). However, when directly injected into an aquifer, nZVI widely interact with bacteria and bacterial biofilms (Goldberg et al., 2007), which may affect nZVI mobility and reactivity. These NP may also present a toxicity for groundwater bacteria. Bacteria present in groundwater are relevant and sensitive indicators of soil or groundwater perturbations because of their key role in biogeochemical cycling (Brookes, 1995; Kandeler et al., 1996; Holden et al., 2014). They can be particularly important in the case of contaminated groundwater because of their ability to metabolize a wide range of contaminants (Antizar-Ladislao and Galil, 2010; Weaver et al., 2015; Zhao et al., 2016). Thus if nZVI application affects bacterial communities then it will also affect biogeochemical reactions occurring in the subsurface.

Several mechanisms of nZVI toxicity towards living organisms have been identified. The first is direct nZVI contact with biological components (O'Carroll et al., 2013). In the majority of studies that used electron microscopy to evaluate damage to cell integrity, precipitation of nZVI or iron oxides on the cell wall or inside the bacterial cells was commonly observed, suggesting that direct contact of nZVI with bacterial cells is required for nZVI to exert toxicity (Auffan et al., 2008; Diao and Yao, 2009; An et al., 2010; Barnes et al., 2010; Li et al., 2010; Xiu et al., 2010; Sevcu et al., 2011; Wille et al., 2017; Jiang et al., 2018; Xue et al., 2018). The second is the release of reactive oxidative species (ROS), such as hydroxyl radicals, superoxide radicals and hydrogen peroxide (Barnes et al., 2010), generated by nZVI in the aqueous phase (Liu et al., 2014). ROS are normally scavenged by antioxidants and various enzymes; however, elevated concentrations of ROS in microbial cells can result in oxidative stress. The damages caused to DNA in the presence of various organisms, including bacteria, were observed in other studies (Sacca et al., 2014; Ghosh et al., 2017). It has also been shown that some bacteria are able to produce proteins that protect them from oxidative stress caused by nZVI (Chaithawiwat et al., 2016; Xue et al., 2018). When nZVI is added into an aqueous medium, a burst of oxidants is produced as Fe0 and ferrous iron (Fe[II]) are converted to ferric iron (Fe[III]). Finally, the release of ferrous iron from nZVI followed by the Fenton reaction can also affect organisms (Xie and Cwiertny, 2012). Furthermore, although the mechanism of nZVI toxicity at the cellular level has not been fully determined, recent studies support the combined effect of biophysical nanoparticle/cell interactions and oxidative stress (Li et al., 2010; Wille et al., 2017; Crampon et al., 2018).

Impact of nZVI towards bacteria has been considerably studied (Jang et al., 2014), sometimes highlighting positive impacts (Wei et al., 2012; Kocur et al., 2016; Dong et al., 2019) and sometimes negative ones (Fajardo et al., 2013; Sacca et al., 2014; Dong et al., 2019). However, these studies generally focus on isolated species rather than addressing the complexity of natural groundwater bacterial communities (Lefevre et al., 2016) although toxicity of xenobiotics may be affected by interactions between the different members of a community. It is thus important to study their impact at the community scale. Indeed, inside a community, the degradation processes are often the result of synergistic degradation or cooperative metabolic interactions (Dejonghe et al., 2003; Freilich et al., 2011), especially in the case of multi-contamination. All members of the community may thus have a key-role in the degradation process. As it has been shown that nZVI may inhibit some groups but not the others inside a bacterial community (Kirschling et al., 2010; Xiu et al., 2010; Fajardo et al., 2013; Xue et al., 2018; Dong et al., 2019), it is important to characterize the toxicity at the scale of the whole community.

In this context, this study aims to deepen our understanding of the impact of nZVI towards groundwater bacterial communities by monitoring the impact of a gradient of nZVI concentrations on the nitrate-reducing potential, the abundance and the structure of a planktonic nitrate-reducing bacterial community. The experimental design was set up in order to monitor bacterial community resistance then resilience after being exposed to nZVI.

### MATERIALS AND METHODS

### Experimental Design

Both the toxicity of nZVI towards bacteria and the resilience of bacteria were studied by using a series of batch experiments. Three distinct steps were carried out and are illustrated in **Figure 1**. Step 1 consisted in enriching a heterotrophic nitratereducing community from a natural groundwater collected from an industrial site by adding a carbon substrate [10 mM sodium acetate (C2H3NaO2)], yeast extracts at 0.2 g L−<sup>1</sup> and sodium nitrate (NaNO3) to a volume of groundwater (v/v) and monitoring nitrate reduction (Wille et al., 2017). Total nitrate reduction was obtained after about 48 h, then experiment proceeded to step 2. Step 2 was carried out for almost 5 days and consisted in two phases. The first one was the inoculation of the pre-culture at a ratio of 1/10 into hermetically sealed 250 mL Schott <sup>R</sup> flasks containing 180 mL of synthetic water enriched with a carbon substrate (10 mM sodium acetate), yeast extracts at 0.2 g L−<sup>1</sup> and sodium nitrate, for a final volume of 200 mL.

The second phase started after 19 h incubation in the dark at 20◦C under agitation and confirmation of triggered nitrate reduction, when nZVI was added to triplicate flasks in order to expose the bacteria to a gradient of nZVI amounts: condition A (control without nZVI), B (8.6 mg Fe.L−<sup>1</sup> ), C (30.8 mg Fe.L−<sup>1</sup> ), and D (70.1 mg Fe.L−<sup>1</sup> ). The addition of nZVI was carried out in an anaerobic chamber under a continuous flow of nitrogen (**Figure 1**). A sterile control without the bacterial inoculum was set up for each condition. 12 mL of liquid samples were collected daily up to 91 h incubation (at 0, 19, 43, 67, and 91 h, samples called T0, T19, T43, T67, and T91, respectively), and analyzed for chemical and biological parameters as described below. nZVI were added to the cultures after the T19 sampling. The third and last step of the experiment consisted of re-inoculating bacterial communities from conditions A, B, C, and D into fresh media (ratio 1/10, synthetic water in the same conditions as step 2) without nZVI to assess community resilience. These new batches were monitored for 3 days until all nitrate was reduced (samples collected at 0, 15, and 39 h incubation).

The synthetic water used in these experiments is FIm water corresponding to moderately hard water as described in US EPA Report EPA-821-R-02-012, Section <IP_ADDRESS> (US EPA, 2002). This water was prepared as previously described in Crampon et al. (2018), and was composed of 96 mg L−<sup>1</sup> NaHCO3, 60 mg L−<sup>1</sup> CaSO4.2H2O, 60 mg L−<sup>1</sup> MgSO4, and 4 mg L−<sup>1</sup> KCl. The pH ranged from 7.4 to 7.8, the hardness ranged from 80 to 100 mg CaCO3.L−<sup>1</sup> and the alkalinity from 57 to 64 mg CaCO3.L−<sup>1</sup> as described by US EPA. The synthetic water was purged with N<sup>2</sup> to remove any oxygen prior to culture medium preparation.

In this experimental setup, as a model activity for the groundwater bacterial community, we focused on bacterial nitrate reduction. This was assessed by measuring nitrate depletion and the abundance of the narG gene and its expression. The overall bacterial community was also characterized by determining abundance (qPCR 16S rRNA), structure (CE-SSCP), composition (Cloning/Sequencing) and viability (LIVE&DEAD <sup>R</sup> ). Physico-chemical parameters were also followed such as pH, redox potential (ORP) and NP size. All analytical approaches are described in the following sections.

### Physico-Chemical Analyses

Fe concentrations were measured by ICP-AES (Jobin Yvon Horiba Ultima 2; Horiba, Japan) in water samples from the cultures at the beginning of the experiment, after acidification with HCl (to obtain a pH < 2). Presented values correspond to the means of two culture replicates for each range of nZVI concentration. Obtained Fe concentrations were 8.6 ± 1.4 mg Fe.L−<sup>1</sup> (n = 2) for the lower concentrations (condition B), 30.8 ± 10.2 mg Fe.L−<sup>1</sup> for intermediate concentrations (condition C) and 70.1 ± 14.8 mg Fe.L−<sup>1</sup> for higher concentrations (condition D).

Due to the presence of ions in the culture medium, NP aggregates were expected to form. The size and zeta potentials of these aggregates were consequently measured. The size of NP aggregates was measured using Dynamic Light Scattering (DLS, Malvern Zetasizer Nano ZS, Malvern, United Kingdom) in cultures at the end of step 2. It was not possible to distinguish here homo- (formed by nZVI aggregation) from hetero- aggregates (formed by aggregation of nZVI with bacteria and/or biofilm). Electrophoretic mobility (EPM) was measured and the apparent zeta potential at 25◦C was calculated automatically by Zetasizer software using the Smoluchowski equation (Hunter, 2001). The polydispersity index (PdI) was used as an indicator of the distribution width. A lower PdI corresponds to a lower sorting of NP, due to the presence of a higher amount of aggregates in solution.

Measurements of pH and redox potential (ORP) were performed for each sampling point in an anaerobic chamber under N<sup>2</sup> atmosphere. Measurements were done with a VTW Multi parameter unit (VTW/Xylem, Inc., Germany) with specific pH and ORP electrodes.

### Nitrate Reduction Activity Measurements

The determination of nitrate concentrations was performed on filtrated water (10 mL, filtration through a 0.22 µm mixed cellulose esters membrane GSWP, Merck Millipore). The analyses were performed at T0, and after 19, 43 h (i.e., 24 h after nZVI addition), 67 and 91 h growth during step 2 and at T0 and after 15 and 39 h growth during step 3. The concentrations of NO3-N were determined by spectrophotometry using a Nova 60 Spectroquant spectrophotometer (Merck/Millipore, Germany) with a specific Spectroquant nitrate kit (reference 114563, Merck/Millipore). The range of concentrations of the kits was 0.1 to 25 mg L−<sup>1</sup> of NO3-N. We decided to focus here on nitrate measurements rather than on the reduction product nitrite as it was only transient and most often below limits of quantification when measured.

# Characterization of Bacterial Communities

### LIVE/DEAD <sup>R</sup> BacLightTM Viability Assays

LIVE/DEAD <sup>R</sup> BacLightTM Bacterial Viability Kits were used according to the manufacturer's instructions (Invitrogen, Carlsbad, CA, United States). In brief, a volume of reactive containing a mixture of SYTO 9 dye (green-fluorescent nucleic acid stain) and Propidium iodide (red-fluorescent nucleic acid stain) is added to an equal volume of culture medium, mixed then incubated in the dark at room temperature for 15 min. Cells are then observed and counted using a fluorescent microscope with filters in the range of 480/500 nm for SYTO 9 stain and 490/635 nm for propidium iodide. SYTO 9 dye tends to stain all bacteria in a population whereas propidium iodide only penetrates cell walls of damaged bacteria, thus showing "live" bacteria as green and "dead" bacteria as red. Observations were made on samples from step 2, at T43, i.e., 24 h after nZVI addition. Preliminary tests (data not shown) were performed on the bacterial community at nZVI concentrations of 10, 30, and 50 mg Fe.L−<sup>1</sup> to perform counts of living and dead cells.

### DNA and RNA Extractions and cDNA Generation

Extractions of DNA and RNA were performed on samples before nZVI addition (T19), after nZVI addition (T43, i.e., 24 h after nZVI addition) and after 67 h incubation (T67, i.e., 48 h after nZVI addition). Biomass was recovered by filtering 10 mL culture on 0.22 µm filters (mixed cellulose esters membrane GSWP, Ø 25 mm, Merck Millipore).

Genomic DNA and total RNA were co-extracted directly from the filters with the combination of MO Bio/QIAGEN (Valencia, CA, United States) kits. First, RNA was extracted with Mo Bio RNeasy <sup>R</sup> PowerSoil <sup>R</sup> Total RNA Kit, and DNA was then recovered thanks to the DNA Elution Accessory Kit, according to the manufacturer's instructions. Remaining DNA was removed from RNA extracts with the turbo DNase (Ambion, Carlsbad, CA, United States) according to the manufacturer's instructions, and RNA was resubmitted to the purification step of the RNeasy <sup>R</sup> PowerSoil <sup>R</sup> Total RNA Kit. The absence of DNA in RNA extracts was confirmed by no PCR signal (see below for PCR conditions). Quantifications of extracted RNA and DNA were performed with the Quantus Fluorometer, using 1 µL extract and the QuantiFluor dsDNA solution (Promega, Charbonnières-les-Bains, France).

The cDNA were obtained from 50 ng RNA extracts by RT-PCR performed on with iScriptTM cDNA synthesis kit and iScriptTM Reverse Transcription Supermix kit (Biorad), following the manufacturer's instructions.

DNA and cDNA were stored at −20◦C and RNA at −80◦C.

### qPCR of 16S rRNA and narG Genes

Quantification of the bacterial 16S rRNA and narG gene copies was performed by qPCR using a CFX ConnectTM Real-Time PCR Detection System (Bio-Rad, France) with primers 341F (3<sup>0</sup> -CCTACGGGAGGCAGCAG-5<sup>0</sup> ) and 515R (30 -ATTACCGCGGCTGCTGGCA-5<sup>0</sup> ) (López-Gutiérrez et al., 2004) for 16S rRNA gene, and primers narG-F (50 -TCGCCSATYCCGGCSATGTC-3<sup>0</sup> ) and narG-R (5<sup>0</sup> - GAGTTGTACCAGTCRGCSGAYTCSG-3<sup>0</sup> ) for nitrate reductase narG gene carried by Proteobacteria (Bru et al., 2007), as only Proteobacteria developed in the batch experiments. qPCR reactions were performed in a total volume of 20 µL, with a master mix containing 7.6 µL of sterile, nuclease- and nucleic acids-free water (MP Biomedicals, Santa Ana, CA, United States), 10 µL of SSO Advanced Universal SYBR Green Supermix (Bio-Rad), 500 nM of each primer, 100 ng of T4 GP32, and 2 µL of template DNA (0.1 to 5 ng µL −1 ) or 2 µL of cDNA. Sterile, nuclease- and nucleic acids-free water was added instead of DNA in no template controls. PCR reactions were performed as follows: 3 min 95◦C, followed by 35 cycles: 30 s 95◦C/30 s 60◦C/30 s 72◦C/for 16S rRNA gene or 6 cycles: 30 s 95◦C/30 s 63◦C/30 s 72◦C, 34 cycles: 30 s 95◦C/30 s 58◦C/30 s 72◦C for narG gene, followed by 30 s at 80◦C for data acquisition, with an additional step rising from 60 to 95◦C at 0.5◦C/s for melt curves generation. All samples, controls and standards were analyzed in duplicates. A calibration curve was obtained from serial dilutions of a known quantity of linearized plasmids containing known copy numbers of 16S rRNA or of narG gene. Results were reported as gene copies per mL of culture medium. Generation of a specific PCR product was confirmed by melting curve analyses and agarose gel electrophoresis.

### CE-SSCP Fingerprints

The bacterial community structure in cultures during the experiment was determined by 16S rRNA gene CE-SSCP fingerprints. About 200 bp of the V3 region of bacterial 16S rRNA genes were amplified from DNA and cDNA with the forward primer w49 (5<sup>0</sup> -ACGGTCCAGACTCCTACGGG-3 0 ; Escherichia coli position, 331) and the reverse primer w34 (50 -TTACCGCGGCTGCTGGCAC-3<sup>0</sup> ; E. coli position, 533), as previously described in Crampon et al. (2018). The 5<sup>0</sup> end was then labeled with the fluorescent dye FAM. 25 cycles and hybridization at 61◦C were used. One µL of diluted (150-fold in nuclease-free water) PCR product was added to a mixture

of 18.8 µL of deionized formamide and 0.2 µL of Genescan-600 LIZ internal standard (Life Technologies, Carlsbad, CA, United States). To obtain single-strand DNA, samples were heat denaturated for 5 min at 95◦C, and immediately cooled on ice. CE-SSCP analyses were performed on an ABI Prism 310 genetic analyzer using a 47-cm length capillary, a nondenaturating 5.6% CAP polymer (Life Technologies) and the following electrophoresis conditions: run temperature 32◦C, sample injection for 5 s at 15 kV, data collection for 35 min at 12 kV. The dataset was analyzed using the software BioNumerics V7.5 (Applied Maths) and lining fingerprints up to the internal standard and to a common baseline.

### Bacteria Taxonomic Identification

About 1,400 bp of the bacterial 16S rRNA gene were amplified from the DNA and cDNA. Universal bacterial primers 8F (5<sup>0</sup> -AGAGTTTGATCMTGGCTCAG-3<sup>0</sup> ) and 1406R (50 -ACGGGCGGTGTGTRC-3<sup>0</sup> ) were used with hybridization at 55◦C and 30 cycles. Purified PCR products (NucleoSpin <sup>R</sup> Gel and PCR clean-up, Macherey–Nagel) were used to construct a gene library with the TOPO-TA Cloning Kit and TOP 10 chemically competent cells (Invitrogen, Carlsbad, CA, United States). To ensure good characterization of the dominant strains detected by CE-SSCP, positive transformants were analyzed by CE-SSCP to find correspondence of clones over the whole community. Based on their migration profile, seven clones were selected for sequencing the insert with vector primers T7 and T3 (Sanger sequencing by GATC Biotech). All seven sequences were submitted to NCBI GenBank under Accession Nos. MH686140 to MH686146.

### Phylogenetic Analysis

The phylogenetic assignment of unambiguous consensus sequences was based on comparison with reference sequences by the NCBI Blast tool<sup>1</sup> . Multiple alignment of sequences and consensus were performed using BioEdit software version <IP_ADDRESS> (Hall, 1999). The phylogenetic tree was constructed using the neighbor-joining method of Jukes and Cantor (Saitou and Nei, 1987), and checked by UPGMA method included in MEGA 10 package (version 10.0.1). The tree topology was tested by bootstrap analysis of 1,000 resamplings. An E. coli 16S rRNA gene sequence was used as an outgroup.

### Statistical Analyses

Statistical analyses were performed on XLStat software version 2017 (Addinsoft) for MS excel. Principal Coordinate Analyses were built from Pearson's correlations coefficients similarity matrices, from triplicates values of the different analyzed parameters. As trends are not linear with nZVI concentrations (inversed trends), control samples were excluded from PCoA analyses.

### RESULTS

### Physico-Chemical Parameters

The mean sizes of nZVI were measured at the end of the first week of the experiment, as well as ZP, in triplicates for each nZVI concentration. ZP of nZVI in cultures was logically negative, and ranged from −37.1 to −23.3 mV.

Globally, PdI was lower for higher nZVI concentrations [PdI = 0.45 ± 0.03 (n = 3) and 0.46 ± 0.08 (n = 3) for Fe concentrations of 30.8 and 70.1 mg Fe.L−<sup>1</sup> , respectively], meaning that a higher concentration of nZVI leads to a lower sorting of aggregates. PdI mean was 0.57 ± 0.11 (n = 3) for lower nZVI concentrations (i.e., 8.6 mg Fe.L−<sup>1</sup> ), showing the presence of less aggregates of different sizes in cultures, and consequently the presence of a higher number of individualized NP.

Even if a trend is observed, showing a higher mean diameter for higher nZVI concentrations (**Table 1**), the mean sizes were not significantly different between the different concentrations (P > 0.05, Kruskal–Wallis test with Conover– Iman pair comparison). Concerning nZVI size sorting shown by PdI values (a value of 1 corresponding to a monodisperse sample, 0 to a polydisperse sample), no significant difference could be highlighted either between the different concentrations (P > 0.05, Kruskal–Wallis test with Conover–Iman pair comparison). However, a trend was again observed (**Table 1**), with a better sorting for lower nZVI concentrations.

Redox potential (ORP) and pH trends were measured in cultures at each sampling date (**Figure 2**). Globally, pH tended to increase during the experiment in all conditions due to the bacterial activity, whereas it was stable in sterile controls (**Figure 2A**). The pH of cultures rose up to around 10 at the end of the experiment. Concerning ORP, it greatly decreased, straight after nZVI addition, logically due to nZVI reduction properties (**Figure 2B**). The lowest ORP values were observed for the highest nZVI concentrations. No ORP decrease was observed just after nZVI addition for the lower concentration (B, 8.6 mg Fe.L−<sup>1</sup> ) compared to the absence of nZVI and to higher nZVI concentrations, but ORP tended to decrease to about −300 mV after 91 h of incubation.

<sup>1</sup>http://www.ncbi.nlm.nih.gov/blast/Blast.cgi

TABLE 1 | Values of Fe concentration, mean NP sizes, PdI and zeta potentials obtained by ICP-AES and DLS analyses in the different cultures.

NA, not applicable.

### Effects of nZVI on Microbial Abundance Molecular Biomass

The total genomic DNA extracted from samples gives an indication of the total molecular biomass. Molecular biomass decreased after nZVI addition for all nZVI concentrations (T43; **Figure 3A**). The decrease was significant (P < 0.05, Student's t-test) for conditions B (8.6 mg Fe.L−<sup>1</sup> ) and C (30.8 mg Fe.L−<sup>1</sup> ), but not for the D modality (70.1 mg Fe.L−<sup>1</sup> ). Moreover, the decrease in molecular biomass was higher for condition B compared to condition C. These results suggest that the highest impact on molecular biomass occurs for lower nZVI concentrations. Forty-eight hours after exposure (T67), molecular biomass increased and reached control values, suggesting a resilience of microorganisms in the community. Surprisingly, at T67, the highest molecular biomass was measured in condition B, where the initial impact on molecular biomass, measured at T43, was highest.

### Bacterial Abundance

As for molecular biomass, nZVI had a significant impact on bacterial abundance, as shown by the decrease in the number of copies of the 16S rRNA gene at T43 in cultures exposed to nZVI in conditions B, C, and D (**Figure 3B**). This impact was significant (P < 0.05, Student's t-test) for all modalities. However, the highest impact of nZVI was again observed for lower concentrations (i.e., 8.6 mg Fe.L−<sup>1</sup> ). Finally, 48 h after exposure (T67), the abundance of the 16S rRNA gene was comparable to values measured in control samples, again suggesting a resilience of the bacterial community.

### Nitrate-Reducing Bacteria Abundance

The determination of narG gene copies gives an indication of the abundance of nitrate-reducing bacteria in the bacterial community. As was observed for molecular biomass (genomic DNA concentration) and bacterial abundance (numbers of 16S rRNA gene copies), a decrease in the abundance of nitrate-reducing bacteria occurred after nZVI addition (T43) independently from nZVI concentrations (**Figure 3C**). This decrease was only significant for conditions B and C (8.6 and 30.8 mg Fe.L−<sup>1</sup> , respectively, P < 0.05, Student's t-test). After 48 h of exposure (T67), the abundance of nitrate-reducing bacteria was comparable to those of the controls. Again, the lowest concentration of nZVI showed the highest impact, but the nitrate-reducing bacteria were resilient, in terms of biomass, after 48 h exposure.

Calculating the ratio between narG and 16S rRNA genes gives an indication of nZVI impact specifically on nitrate-reducing bacteria over the whole community (**Figure 3D**). Before nZVI addition, about 90% of total bacteria carried the narG gene. This value is significantly lower after nZVI addition for conditions B and C, especially for the lowest nZVI concentration. This value drops to around 70% of nitrate-reducing bacteria in the whole community for condition B, and 80% for condition C. The ratio of nitrate-reducing bacteria is the same for condition D (70.1 mg Fe.L−<sup>1</sup> ) compared to control. Overall, these parameters remained stable for the control condition A.

### Bacteria Viability

In the control condition A (no nZVI), many live bacteria aggregates were observed after 43 h of incubation [T43; A(1), A(2), and A(3) in **Figures 4a,A**]. However, with nZVI (**Figures 4a,B–D**), bacterial aggregates were not found, except smaller ones in condition D, and more dead cells were present.

A viability of around 80% in the control condition was observed (counts were performed on diluted samples with between 10 and 50 cells per image). This value drops to around 5% in condition B, 20% in condition C and 40% for the highest concentration in condition D (**Figure 4b**). The higher toxicity in the presence of lower concentrations of nZVI is thus again confirmed here by LIVE/DEAD <sup>R</sup> observations, with a smaller number of living cells in the case of lower nZVI concentrations.

FIGURE 4 | LIVE/DEAD <sup>R</sup> observations (illustration images) (a) and results of counting (b) of bacteria after 43 h of incubation (24 h after nZVI addition).  Please note that images from (a) were not those used for LIVE/DEAD counting results presented in (b).

# Effects of nZVI on Microbial Activity

### Nitrate Reduction Kinetics

fmicb-10-00533 March 19, 2019 Time: 11:44 # 8

The nitrate concentrations in sterile controls remained stable all along the experiment. Moreover, nZVI had no significant impact on nitrate concentrations (**Supplementary Material SM2**). Therefore, the observed nitrate reduction in cultures is of bacterial origin. Before nZVI addition (**Figure 5A**), nitrate reduction occurred and was repeatable for all conditions. After nZVI addition, nitrate reduction stopped in all cultures, except in control without nZVI. However, 24 h after nZVI addition (T43), nitrate reduction activity was again observed in all cultures with nZVI. After 91 h, remaining nitrate concentrations were even lower than in cultures without nZVI. These results highlight a resilience capability of the nitrate-reducing bacterial community after exposure to nZVI, even if nZVI initially had a significantly toxic impact on nitrate reduction (P > 0.05, Kruskal–Wallis test with Conover–Iman pair comparison).

During the resilience step (second week), nitrate reduction was monitored until complete nitrate disappearance (**Figure 5B**). For all cultures previously exposed to nZVI during the first week experiment, a 15-h lag-phase was observed. The difference was significant between cultures that had not been exposed to nZVI (A), for which no lag phase was observed, and cultures that had been exposed to nZVI (B, C, and D), independently from nZVI concentrations (P > 0.05, Kruskal–Wallis test).

Of note, in complement to sterile controls presented here, a preliminary experiment was led, with nZVI concentrations ranging from 1 to 500 mg Fe.L−<sup>1</sup> in sterile conditions and with comparable NO<sup>3</sup> concentrations, to study the impact of nZVI on nitrates (**Supplementary Material SM1**). The presence of nZVI did not impact nitrate concentrations, independently from nZVI concentrations.

### Bacterial Activity: RNA vs. DNA

As ribosome content (rRNA) in bacterial cells is an indicator of active or recently active bacteria, ratios of 16S rRNA copy numbers (determined from cDNA generated from total RNA) to 16S rRNA gene copy numbers (determined from DNA) was used to evaluate the growth state of the bacterial community (**Figure 6A**). 16S rRNA/16S rDNA ratio was highest 24 h after nZVI addition (T43) compared to before addition (T19) in all conditions. It was significantly higher (P-value < 0.05) for conditions B and D, and, although statistically un-significant, the same trend was observed for conditions A and C. Thus, despite a negative effect on microbial biomass and bacterial abundance, nZVI did not affect the global growth state of bacteria remaining in the community.

### Nitrate Reduction Activity: narG mRNA/16S rRNA and narG mRNA/narG DNA

narG transcripts (narG mRNA) were used to evaluate the nitrate reduction activity. The ratio between the number of narG mRNA and 16S rRNA (**Figure 6B**) highlights the part of nitrate-reducing activity among active or recently active bacterial cells in the community. No significant changes were observed between controls and cultures containing nZVI or between samples collected just before (T19) and after (T43) nZVI addition. Therefore, the impact of nZVI on the nitrate reduction activity was not clearly highlighted here. Only trends were observed at T43, with a decrease in narG transcripts/16S rRNA ratios, suggesting a decrease in the part of nitrate reduction among active or recently active bacterial cells, and again for lower nZVI concentrations. The ratio of narG mRNA to narG DNA shows the nitrate reduction activity of the nitrate reduction community. No significant differences were observed between controls and conditions with nZVI, or between samples collected just before (T19) and 24 h after (T43) nZVI addition, suggesting that the remaining nitrate-reducing bacteria were active. Even if not significant, the narG mRNA/DNA ratio tends to increase at T43 for cultures with nZVI, suggesting that nZVI addition, especially at lower concentrations, slightly increased the activity of remaining nitrate-reducing cells (**Figure 6C**).

In order to compare the impact of the different nZVI concentrations on biological and physico-chemical parameters at the different times, statistical analyses were performed (PCoA, **Supplementary Material SM3**). No correlations between nZVI concentrations and the different studied parameters were observed, as expected, before nZVI addition (T19). After nZVI addition, the maximum correlation was observed at T43 (24 h after nZVI addition), but was no longer relevant at T67, indicating that the toxicity of nZVI, for these parameters, lasted for no more than 48 h.

# Effects of nZVI on Microbial Community Structure

Beyond the impact on abundance and activity, the bacterial community structure could also be affected by the presence of nZVI. Community CE-SSCP fingerprints of the 16S rRNA gene were performed on (i) DNA and (ii) cDNA, which allows to observe (i) the structure of the total community and (ii) the structure of active or recently active communities. Before nZVI addition, fingerprints showed three main peaks with one dominating the others, and suggested a low bacterial diversity in all conditions. Overall, fingerprints were similar for DNA and cDNA, highlighting a good adequacy between detection and growth state of all detected bacteria in the cultures. After nZVI addition, a clear shift of the bacterial community structure occurred. A new peak (in red) was detected to the detriment of the dominant peak (in green in **Figure 7**), which was highly affected by nZVI addition (**Figures 7B–D**).

Indeed, especially at T67, a shift in diversity was observed, with the detection of a new peak (in red in **Figure 7**) at the expense of the main peak (in green). This shift was more marked for lower nZVI concentrations (**Figure 7B**), where the new peak became dominant on the fingerprints (**Figure 7**). However, as nZVI concentrations increase, the new peak appears less and less dominant, and proportions between green and red peaks are quite similar for conditions D (**Figure 7**). After nZVI addition, the relative proportion of bacteria corresponding to peaks 1 and 2 increased in the community. This result is intermediate for condition C (**Figure 7C**). For condition B, the new peak in red

appeared at T43, 24 h after nZVI addition, and became dominant at T67, 48 h after nZVI addition. Finally, the relative proportion of increasing bacteria, corresponding to the red peak, was the highest at T67, especially for B and C conditions for which it became the dominant peak. Globally, nZVI had an impact on the bacterial community structure, either favoring or affecting presence and abundance of the bacterial strains.

### Bacteria Identification

Cloning and sequencing of community 16S rRNA gene at the end of the experiment (T67) and CE-SSCP migration patterns of retrieved clones allowed to attribute several sequences to the main peaks of the fingerprints (**Figure 8A**). All 16S rRNA gene clone sequences were closely related (>99% blast identity) to that of members of the Pseudomonas genus (**Figure 8B**). Species from the Pseudomonas genus are difficult to distinguish with 16S rRNA gene sequences, as they are very closely related to each other (Anzai et al., 2000). However, although we cannot clearly link a sequence to a species, phylogenetic comparison with close members of the Pseudomonas genus confirmed the shift in community composition. For the first peak of interest (in green), prior to nZVI addition, three clones were identified (Clone\_nZVI\_Ps#1, #2, and #3). In particular, Clone\_nZVI\_Ps#1 branched among facultative anaerobic Pseudomonas that are able to perform nitrate reduction (**Supplementary Material SM4**). Consequently, Pseudomonas strains corresponding to the dominant green peak in the fingerprints could be the one able to perform nitrate reduction, and could be involved in the nitrate reduction process observed in the experiments. For the second peak of interest (in red), after nZVI addition, the three clones identified (Clone\_nZVI\_Ps#4, #5, and #6) are phylogenetically different from Clone\_nZVI\_Ps#1 and #3. They are highly related to each other (dissimilarity < 0.01%), and very close to Pseudomonas kunmingensis, Pseudomonas knackmussii, and Pseudomonas chloritidismutans sequences. Pseudomonas kunmingensis is an aerobic bacteria able to perform nitrate reduction (Romanenko et al., 2005; Xie et al., 2014). Pseudomonas knackmussii and Pseudomonas chloritidismutans are facultative aerobic bacteria, unable to perform nitrate reduction (**Supplementary Material SM4**). However, as nitrate reduction then resumed, either the rising strain could reduce nitrate or nitrate reduction involved bacteria not detected by CE-SSCP. Indeed, one last clone sequenced (Clone\_nZVI\_Ps#7) also affiliated to Pseudomonas, did not correspond to a peak visible on the CE-SSCP fingerprint.

### DISCUSSION

This work studied the impact(s) of nZVI addition on the abundance, the activity and the structure of a bacterial community grown from a natural groundwater. This was experimentally assessed by focusing on a key function in groundwater, nitrate reduction. Moreover, a follow-up experiment was carried out to look at the bacterial community's reaction when inoculated in fresh culture medium after their exposure to nZVI.

Addition of nZVI, even at low concentrations of ∼9 mg L−<sup>1</sup> significantly decreased molecular biomass and bacterial abundance 24 h after addition as shown by LIVE/DEAD <sup>R</sup> observations, molecular biomass and 16S rRNA

and narG gene abundances. However, molecular biomass then increased again rapidly in all conditions suggesting only a temporary effect on bacterial abundance. Few studies have assessed the impact of nZVI in planktonic groundwater or water studies. Barnes et al. (2010) observed an increase in CFU counts in river water submitted to 100 mg L−<sup>1</sup> nZVI with no initial decrease. Kirschling et al. (2010) also observed no impact of 1.5 g L−<sup>1</sup> nZVI on the number of 16S rRNA gene copies in bacterial communities from TCE contaminated groundwater. However, this was shown after 250 days incubation so it is possible that the community had time to adapt. Fajardo et al. (2012) found little impact of 34 mg.nZVI.g−<sup>1</sup> soil on microbial

viability; however, as underlined by Fajardo et al. (2013) and Sacca et al. (2014), the impact of nZVI in soils is closely linked to their composition whereas in the present study, in an aqueous solution, bacterial communities would be more directly exposed, as it would be the case in aquifers (less complex compared to soils, low content of natural organic matter).

Similarly, the growth state and nitrate-reducing activity, assessed by the ratio of narG transcripts and nitrate depletion, were only affected temporally by the three concentrations of nZVI, suggesting an adaptation of the bacterial community to the new environment containing nZVI. The activity of remaining bacteria in condition B, the most affected by nZVI addition,

even appeared to be stronger compared to the less impacted cultures. At the end of the experiment, no nitrates were left in nZVI-exposed cultures, whereas a residue was always found in controls without nZVI. One explanation could be anaerobic corrosion of nZVI, leading to the release of H<sup>2</sup> in the media, a well-known electron donor for bacteria, reinforcing their metabolism and resulting, after an adaptation period, to a more complete nitrate reduction in the presence of nZVI. A lower nitrate concentration was also observed during the resilience step for cultures having been exposed to nZVI, and could be explained both by residual H<sup>2</sup> introduced in fresh media by the inoculum (1/10 volume), or by the bacterial community structure modifications generated by nZVI during the first step of the experiment. Even if nitrate reduction activity resumed and if parameters related to bacterial abundance and activity returned to "normal" values, the structure of the bacterial community was modified by nZVI. Indeed CE-SSCP fingerprints revealed the increase of a new peak after nZVI addition, first in condition B which was the most impacted globally, and then in conditions C and D. This new peak then dominated the incubations and was confirmed to be a different strain to the initially dominating peak by sequencing fragments of 16S rRNA gene. These changes in the community structure due to nZVI were previously studied, and attributed both to H<sup>2</sup> release and changes in redox conditions, favoring the growth of methanogen or sulfate-reducing bacteria (Kirschling et al., 2010; Xiu et al., 2010), but also to some specific nZVI coatings, protecting bacteria from oxidative stress, and creating a new carbon source in the media (Kirschling et al., 2010). In our case (non-coated nZVI), a rapid impact could be observed, attributable to toxic effect (oxidative stress), but the community structure was definitely modified even if nitrate-reduction resumed after a lag-phase This impact on community structure is highly dependant on the media, as was shown in soils (Fajardo et al., 2018). The variability of the media is however lower in groundwater compared to soils, and so changes in community structure due to nZVI are expected to occur more currently. These results are finally comparable to several previous experiments carried out either in water samples or soils. Indeed, both Fajardo et al. (2012) and Sacca et al. (2014) observed significant changes in bacterial community diversity and structure when applying nZVI to soils as well as Kirschling et al. (2010) in contaminated groundwater using a DGGE approach, whereas Barnes et al. (2010) found no impact of nZVI in river-water bacterial communities.

When re-cultivated in fresh media without nZVI (step 3), the temporary inhibition in nitrate reduction observed for cultures previously exposed to nZVI could also correspond to an adaptation period in the opposite way, due to the absence of nZVI in fresh substrate. Lag phase represents the earliest and most poorly understood stage of the bacterial growth cycle (Rolfe et al., 2012). In bacterial growth models, lag phase is generally described as an adaptation phase of the bacterial community when introduced in a new environment. In our case, the changing parameter is the absence of nZVI in the fresh culture media (lower oxidative stress, lower Fe concentrations). Bacterial growth after this lag phase is expected to be exponential (Rolfe et al., 2012). As the conditions were the same for the cultures without nZVI, it therefore seems logical that no lag phase was observed. As previously shown (Rolfe et al., 2012), high concentrations of iron could be observed in bacterial cells during lag phase, and these concentrations were associated with transient sensitivity to oxidative stress. Moreover, it was shown that toxicity of nZVI is growth phase dependent (Chaithawiwat et al., 2016). A new exposition to nZVI during this lag phase (corresponding to a resilience period for the bacterial communities), would possibly have been lethal for the bacterial community.

Surprisingly, toxicity was higher for lower nZVI concentrations. Even if DLS analyses did not enable to highlight significant differences between the different concentrations in terms of mean sizes of NP aggregates and NP size sorting, trends were observed: more single nZVI or small aggregates were present in lower nZVI concentrations, and higher nZVI concentrations tend to result in the formation of higher aggregates. Consequently, the increase of NP concentrations enhances their aggregation. These observations are consistent with the literature (Baalousha, 2009). During exposure to nZVI, pH tended to increase, and was always above pH 8.0. For pH > 7, functional groups of bacteria and biofilms are globally deprotonated, and are thus negatively charged (Fein et al., 2001; Borrok and Fein, 2004). The formation of ionic bridges between negatively charged bacterial cells and NP with divalent cations present in solution (e.g., Ca2<sup>+</sup> and Mg2+), could lead to strong interactions between bacterial membranes and nZVI, and thus to toxic effects. The presence of large aggregates could result in a decrease in these interactions, and consequently in the observed decrease in toxicity. The formation of large aggregates in cultures with higher nZVI concentrations decreases the ZVI area per unit volume of liquid (Baalousha, 2009), and could also lead to the observed lower toxicity for higher nZVI concentrations. ROS generated would indeed be lower in the media with lower active ZVI area per volume unit, leading to a lower toxic effect.

As toxicity may be due to different mechanisms, these mechanisms may differ as a function of nZVI concentrations. At high concentrations, the main toxicity mechanism could be due to the release of ROS. At low concentrations, it could be rather due to cell/NP interactions. There is probably a cooccurrence of these mechanisms in any case, but their relative part could be dependent on nZVI concentrations and subsequent aggregation. As Fajardo et al. (2012) suggested for soils, these results point to a dose and species-specific impact of nZVI applications in groundwater.

### CONCLUSION

This study aimed to characterize the impact of nZVI addition on a bacterial community grown in controlled conditions from a natural groundwater sample. The impact of nZVI addition was assessed on bacterial abundance, activity (especially denitrifying activity) and community structure. Even if the impact is only temporary on abundance and activity, due to a good resilience/adaptation of the bacterial community, a shift was observed in community structure. Remediation of contaminated

groundwater using nZVI could consequently expose the bacterial groundwater community to shifts in their structures, and thus such techniques could have a long-term impact on bacterial communities and by extension on groundwater ecological functions, such as nitrate reduction. It is consequently important to study the impact on the bacterial community structure and on potential metabolic functions of the aquifers, at the long-term and at the field-scale.

### AUTHOR CONTRIBUTIONS

MCr contributed to the experimental design, carried out experimental work, data analysis, and drafted the manuscript. JH, CJ, and PO contributed to the experimental design, data analysis, and drafting and revising the manuscript. MCh carried out some of the experiments linked to molecular biology in the study.

### FUNDING

We are thankful for funding received from the European Union Seventh Framework Programme (FP7/2007–2013) under Grant

### REFERENCES

Agreement No. 309517. We gratefully acknowledge the financial support provided to the PIVOTS project by the Région Centre – Val de Loire (ARD 2020 program and CPER 2015–2020) and the French Ministry of Higher Education and Research (CPER 2015–2020 and public service subsidy to BRGM). This operation was co-funded by European Union. Europe was committed to the Centre-Val de Loire region with the European Regional Development Fund.

### ACKNOWLEDGMENTS

We thank NANOIRON (Rajhrad, Czechia) for providing the nZVI used in the present experiments. We also thank the Carnot program for financing the post-doctoral position that allowed conducting this study.

### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2019.00533/full#supplementary-material

phase dependent. Chemosphere 144, 352–359. doi: 10.1016/j.chemosphere.2015. 09.025

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2019 Crampon, Joulian, Ollivier, Charron and Hellal. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Benthic Diatom Communities in an Alpine River Impacted by Waste Water Treatment Effluents as Revealed Using DNA Metabarcoding

Teofana Chonova1,2 \*, Rainer Kurmayer<sup>1</sup> , Frédéric Rimet<sup>2</sup> , Jérôme Labanowski<sup>3</sup> , Valentin Vasselon<sup>2</sup> , François Keck2,4, Paul Illmer<sup>5</sup> and Agnès Bouchez<sup>2</sup>

<sup>1</sup> Research Department for Limnology, Mondsee, Faculty of Biology, University of Innsbruck, Mondsee, Austria, <sup>2</sup> UMR CARRTEL, INRA, Université Savoie Mont Blanc, Thonon-les-Bains, France, <sup>3</sup> UMR IC2MP 7285, CNRS, Université de Poitiers, ENSIP, Poitiers, France, <sup>4</sup> Department of Aquatic Sciences and Assessment, Faculty of Natural Resources and Agricultural Sciences, Swedish University of Agricultural Sciences, Uppsala, Sweden, <sup>5</sup> Department of Microbiology, Faculty of Biology, University of Innsbruck, Innsbruck, Austria

### Edited by:

Ghiglione Jean-Francois, Centre National de la Recherche Scientifique (CNRS), France

### Reviewed by:

Gabor Varbiro, Centre for Ecological Research (MTA), Hungary Yonghong Bi, Institute of Hydrobiology (CAS), China

> \*Correspondence: Teofana Chonova<EMAIL_ADDRESS>
### Specialty section:

This article was submitted to Microbiotechnology, Ecotoxicology and Bioremediation, a section of the journal Frontiers in Microbiology

Received: 15 December 2018 Accepted: 15 March 2019 Published: 09 April 2019

### Citation:

Chonova T, Kurmayer R, Rimet F, Labanowski J, Vasselon V, Keck F, Illmer P and Bouchez A (2019) Benthic Diatom Communities in an Alpine River Impacted by Waste Water Treatment Effluents as Revealed Using DNA Metabarcoding. Front. Microbiol. 10:653. doi: 10.3389/fmicb.2019.00653 Freshwater ecosystems are continuously affected by anthropogenic pressure. One of the main sources of contamination comes from wastewater treatment plant (WWTP) effluents that contain wide range of micro- and macropollutants. Chemical composition, toxicity levels and impact of treated effluents (TEs) on the recipient aquatic ecosystems may strongly differ depending on the wastewater origin. Compared to urban TEs, hospital ones may contain more active pharmaceutical substances. Benthic diatoms are relevant ecological indicators because of their high species and ecological diversity and rapid response to human pressure. They are routinely used for water quality monitoring. However, there is a knowledge gap on diatom communities' development and behavior in treated wastewater in relation to prevailing micro- and macropollutants. In this study, we aim to (1) investigate the response of diatom communities to urban and hospital TEs, and (2) evaluate TEs effect on communities in the recipient river. Environmental biofilms were colonized in TEs and the recipient river up- and downstream from the WWTP output to study benthic diatoms using DNA metabarcoding combined with high-throughput sequencing (HTS). In parallel, concentrations of nutrients, pharmaceuticals and seasonal conditions were recorded. Diatom metabarcoding showed that benthic communities differed strongly in their diversity and structure depending on the habitat. TE sites were generally dominated by few genera with polysaprobic preferences belonging to the motile guild, while river sites favored diverse communities from oligotrophic and oligosaprobic groups. Seasonal changes were visible to lower extent. To categorize parameters important for diatom changes we performed redundancy analysis which suggested that communities within TE sites were associated to higher concentrations of beta-blockers and non-steroidal anti-inflammatory drugs in urban effluents vs. antibiotics and orthophosphate in hospital effluents. Furthermore, indicator species analysis showed that 27% of OTUs detected in river downstream communities were indicator for urban or hospital TE sites and were

absent in the river upstream. Finally, biological diatom index (BDI) calculated to evaluate the ecological status of the recipient river suggested water quality decrease linked to the release of TEs. Thus, in-depth assessment of diatom community composition using DNA metabarcoding is proposed as a promising technique to highlight the disturbing effect of pollutants in Alpine rivers.

Keywords: pharmaceuticals, diatom communities, WWTP effluents, DNA metabarcoding, functional traits, water quality index, indicator species analysis

### INTRODUCTION

Freshwater ecosystems provide important resources and services to humans. However, their sustainability is nowadays constantly affected due to their permanent application for conflicting purposes (e.g., release of wastewater and production of drinking water). WWTPs are used to reduce the release of anthropogenic pollutants into aquatic ecosystems, but pollutants cannot be completely eliminated by the treatment process (e.g., Verlicchi et al., 2015). Thus, TE usually contains a wide spectrum of highly concentrated macro- and micropollutants that may threaten ecosystems health (e.g., Verlicchi et al., 2012; Chonova et al., 2017).

Composition and toxicity level of effluents may differ strongly depending on their origin (e.g., urban, hospital, industrial), and may thereby have different impact on aquatic ecosystems (Labanowski et al., 2016). As recently shown by Chonova et al. (2016), microbial biofilm communities from urban and hospital TE sites exhibit remarkable differences in their development. Thus, the origin of pollution should not be neglected in the context of aquatic environmental monitoring. Increasing anthropogenic pressure on freshwater ecosystems imposes politically defined restrictions linked to their ecological status and functioning. In Europe, the water framework directive (WFD) aims at ensuring the maintenance of good water quality by monitoring chemical hazard substances and key aquatic indicator organisms (European Commission, 2013).

Benthic diatoms are microalgae used worldwide for water quality assessment. They are of particular interest in the context of bioassessment because of their taxonomic diversity and different species sensitivity and resistance to pollution. Diversity and community composition of diatoms adapt rapidly to the presence of chemical, physical, and biological disturbances (e.g., Stevenson and Smol, 2003). In the last decades, these features were used for the development of diatom indices that rely mainly on taxonomic composition (species or genera) obtained from microscopy counts and serve to evaluate general pollution (e.g., Coste et al., 2009; Schneider and Lindstrøm, 2009, 2011). The efficiency of DNA metabarcoding combined with high-throughput sequencing (HTS) techniques in this context was recently largely explored (e.g., Visco et al., 2015; Zimmermann et al., 2015; Vasselon et al., 2017a,b). The ability of these innovative molecular methods to provide fast, costefficient and reliable diatom inventories makes them promising biomonitoring tool and their application is currently being improved (e.g., Vasselon et al., 2017a, 2018).

Several studies have suggested that diatom sensitivity to environmental pressures is closely related to ecological guilds sharing functional traits (Passy, 2007; Tapolczai et al., 2016, 2017) which is also reflected by the phylogenetic tree of diatoms (Keck et al., 2016; Esteves et al., 2017). These relations give complementary insights into the health of aquatic ecosystems (Rimet and Bouchez, 2012; Larras et al., 2014) avoiding drawbacks of precise taxonomic identification (e.g., ecoregional differences in species ecological optima, taxonomic misidentification). However, the taxonomic and ecological complexity requires better understanding and further investigations.

In the context of bioassessment, diatom sensitivity is principally explored and applied in relation to nutrients and organic matter concentrations. However, the high variety of micropollutants constantly released in aquatic environments raise serious concerns about the health of these ecosystems (Schwarzenbach et al., 2006). Besides macropollutants such as nitrogen, phosphorus, and dissolved organic carbon (DOC), the possible interaction and impact of micropollutant mixtures on diatom communities should not be neglected. The sensitivity of diatoms to micropollutants (e.g., pesticides, heavy metals, and pharmaceuticals) at environmentally relevant concentrations was reported in numerous ecotoxicological studies. However, considering the large variety of diatoms and chemicals in the environment, single-species bioassays are not sufficient to provide environmentally realistic picture (Hagenbuch and Pinckney, 2012). Studies at higher complexity level are useful to fill knowledge gaps on toxic impacts. Rimet and Bouchez (2012) for example studied life-forms, ecological guilds, and cell size of diatoms in mesocosm experiments to develop a tool to assess pesticide contamination. Larras et al. (2014) found a relation between diatom sensitivity to herbicides and species' phylogenetic position. However, assessing impacts of micropollutant mixtures on the entire diatom community in natural ecosystems, and distinguishing them from the effects induced by other factors (e.g., nutrients, flow velocity, light) remains challenging (Marcel et al., 2013).

**Abbreviations:** BDI, biological diatom index; COD, chemical oxygen demand; H, hospital treated effluent; HPLC–MS/MS, high performance liquid chromatography coupled to a mass spectrometer; HTS, high-throughput sequencing; ISA, indicator species analysis; NMDS, non-metric multidimensional scaling; NSAID, nonsteroidal anti-inflammatory drug; OTU, operational taxonomic unit; PCR, polymerase chain reaction; PERMANOVA, permutational multivariate analysis of variance; PhC, pharmaceutical compounds; pRDA, partial redundancy analysis; RD, river downstream; RDA, redundancy analysis; RU, river upstream; TSS, total suspended solids; TEs, treated effluents; U, urban treated effluent; VIF, variance inflation factor; WFD, water framework directive; WWTP, wastewater treatment plant.

Tapolczai et al. (2016) suggests that habitats with increased and multiple stress conditions are well adapted for studying the response of benthic diatom communities to environmental pressure. Modifications of diatom communities in rivers downstream from WWTP outputs have been already reported (e.g., Tornés et al., 2018), but community dynamics in direct response to various TEs remained thereby unknown. Despite the regular use of diatoms as biological indicators, there is a knowledge gap on their development and behavior in treated WWTP effluents. Understanding the dynamics of such communities is important to disentangle toxicity effects of pollutants mixtures and to better define a possible role of diatoms as bioindicators. These communities can be seen as a reference to better understand and evaluate the impact of TEs on recipient rivers.

In this study, we investigated community dynamics of periphytic diatoms that were exposed to hospital (H) and urban wastewater treated effluents (U), and in the recipient river up- (RU) and downstream (RD) from the WWTP output. Thereby, we aimed to (1) compare the effect of two highly polluted environments with different concentrations of pharmaceutical compounds (PhC) on diatom communities, and (2) evaluate and compare the effect of TEs' release on benthic diatom community composition in the recipient river. Diatom communities were analyzed using a DNA metabarcoding approach. Characterization of communities regarding their diversity, taxonomic composition, phylogeny, and functional traits were used to give complimentary insights on diatoms development and behavior in habitats with different loads of pollution. Indicator species analysis (ISA) and biological diatom index (BDI) were used to evaluate the influence of TEs' release on river communities.

### MATERIALS AND METHODS

# Study Site and Biofilm Colonization Experiments

The experimental site is a pilot WWTP handling separately urban and hospital wastewater using biological treatment with conventional activated sludge system (Chonova et al., 2018b) (**Figure 1A**). The hospital whose effluents flow directly without special pretreatment to the WWTP has 450 beds, and the urban network includes around 20,850 inhabitants. This parallel separate treatment enables comparison of urban and hospital TEs and their effect on benthic diatoms. Treated wastewaters are discharged into the recipient River Arve. Downstream from the WWTP (ca. 18 km), water from River Arve is used for drinking water production for the city of Geneva (Switzerland), which increases the importance for the monitoring and maintenance of its high water quality.

The survey was performed between February and July 2014 in the following four locations: urban (U) and hospital (H) treated effluents and recipient river approximately 250 m

FIGURE 1 | (A) Sampling site map (U, urban treated effluent; H, hospital treated effluent; RU, river upstream; RD, river downstream; WWTP basins are indicated with light-blue circles; approximate sampling locations are labeled with dots in the respective color); (B) sampling location river; (C) WWTP basin; (D) sampling location treated effluent; (E) blueschist stones exposed for biofilm colonization (approximate stone surface: 80 ± 20 cm<sup>2</sup> ).

upstream (RU) and 40 m downstream (RD) from the WWTP output (**Figure 1**). Differences in hydraulic conditions between river and TE sampling locations were observed, suggesting a higher flow velocity and turbulence in river vs. laminar conditions exhibiting at least twice lower velocity in TE sampling locations. To study the development of benthic diatom communities in each location and compare them to natural seasonal gradient, biofilm colonization experiments were done as described in Chonova et al. (2016). Metal grid-baskets with clean and previously autoclaved blueschist stones (approximate stone surface: 80 ± 20 cm<sup>2</sup> ) were installed as substrates for natural biofilm colonization in water depth between 10 and 40 cm and were regularly inspected. Total surface of stones sampled in each location was 200 cm<sup>2</sup> at least which is in accordance with the French norms used for application of WFD. After each colonization period (about 1 month), the biofilm was scraped from the stones and suspended in sterile water. Biofilm samples were collected in triplicates obtained from independent stones to ensure reliability of the experimental design. Samples were then immediately transported to the laboratory in cooling boxes for further analysis. This colonization experiment was repeated six times between February and July 2014.

## Characterization of Micro- and Macropollutants and Seasonal Conditions

Sampling sites were previously characterized in terms of nutrients (ammonium, nitrite, nitrate, orthophosphate), TSS, COD, and set of pharmaceuticals from the therapeutic classes beta-blockers (atenolol, propranolol), NSAIDs (diclofenac, ibuprofen, and ketoprofen), antibiotics (ciprofloxacin, sulfamethoxazole, and vancomycin), analgesics (paracetamol), and anticonvulsants (carbamazepine) as described in Chonova et al. (2016, 2018a). Flow-proportional 24-h sampling was performed. To obtain a representative flow-proportional sample, first, time proportional subsamples were collected every 10 min during 24 h. Secondly, subsamples from the same hour were pooled together to obtain one time-proportional sample per hour (resulting in 24 time-proportional samples of 1,020 mL in total). Finally, a fraction of each of them was combined according to the water flow rate. Nutrients, TSS and COD were measured with standard analytical methods, described by the French standard operating procedures (AFNOR, 1997). Pharmaceuticals were measured by solid-phase extraction (SPE) using hydrophilic–lipophilic balanced (HLB) columns and analyzed by HPLC–MS/MS as described by Wiest et al. (2018). For each parameter and sampling location, we calculated mean concentrations (from six monthly measurements for TE sites, and from three measurements performed in November and December 2013 and March 2014 for river sites) to give an overview of the different habitat characteristics. Measurements of hydro-meteorological conditions (river flow, urban and hospital WWTP discharge, solar irradiance, river and air temperature and precipitation) were collected as described in Chonova et al. (2018a) and means from multiple daily measurements were calculated for each colonization period to characterize seasonal trends. A theoretical dilution factor for the RD site was calculated by dividing the average river flow (m<sup>3</sup> .d) by the average WWTP discharge (m<sup>3</sup> .d−<sup>1</sup> ).

### Metabarcoding of Diatom Samples DNA Extraction and Polymerase Chain Reaction (PCR) Amplification

Total genomic DNA was isolated using Sigma-Aldrich GenEluteTM-LPA as described in previous studies (e.g., Chonova et al., 2016). This method was recommended for diatom metabarcoding, because of its combination of various lysis mechanisms that are helpful for diatom cells opening and its ability to provide a large quantity of DNA (Vasselon et al., 2017a).

A 312 bp fragment of the rbcL plastid gene, recommended as DNA barcode for diatom metabarcoding (Kermarrec et al., 2013), was amplified using Takara LA Taq <sup>R</sup> polymerase and the primer pairDiat\_rbcL\_708F (AGGTGAA GTTAAAGGTTCATACTTDAA) (Stoof-Leichsenring et al., 2012) and R3 (CCTTCTAATTTACCAACAACTG) (Bruder and Medlin, 2007). Each PCR amplification mix (total volume of 25 µL) contained 1 µL of 25 ng DNA template, 0.75 U of Takara LA Taq <sup>R</sup> polymerase, 2.5 µL of 10× PCR Buffer, 1.25 µL of 10 µM of each primer, 1.25 µL of 10 g/L BSA, 2 µL of 2.5 mM dNTP, and ultrapure water to complete. For each set of reactions, a negative control was included. PCR reaction conditions were initiated by a denaturation step at 95◦C for 15 min followed by a total of 30 cycles of 95◦C for 45 s (denaturation), 55◦C for 45 s (annealing), and 72◦C for 45 s (final extension) (Vasselon et al., 2017a).

### Genomic Libraries Preparation

Genomic libraries for HTS were prepared from PCR amplicons of the rbcL 312 bp fragment from four independent PCR reactions per sample as described in Vasselon et al. (2017a). Briefly, PCR amplicons were cleaned with Agencourt AMPure beads (Beckman Coulter, Brea, CA, United States). Quality and quantity of purified products were measured with 2200 Tape Station (Agilent Technologies, Santa Clara, CA, United States). Tags were added to each amplicon and libraries were prepared using the NEBNext <sup>R</sup> Fast DNA Library Prep set for Ion TorrentTM (BioLabs, Ipswich, MA, United States) and A-X tag adapter provided by Ion ExpressTM Barcode adapters (Life Technologies, Carlsbad, CA, United States). Finally, all libraries were pooled at a final concentration of 100 pM and sequenced using an Ion 316TM Chip Kit V2 (Life Technologies, Carlsbad, CA, United States) on a PGM Ion Torrent sequencer by the "Plateforme Genome Transcriptome" (PGTB, Bordeaux, France).

### Sequence Data Processing

The sequencing platform performed demultiplexing and adapter removal. Data were provided in form of separate fastq files for each sample. Bioinformatic processing was performed with Mothur Software (Schloss et al., 2009) as proposed by Vasselon et al. (2017a). Low quality sequences (read length below 250 bp, Phred quality score below 23 over a moving window of 25 bp,

more than one mismatch in the primer sequence, homopolymer over 8 bp or presence of ambiguous base) were removed and obtained quality filtered sequences from all samples were analyzed jointly. Removal of chimera using Uchime algorithm (Edgar et al., 2011) was performed after pre-clustering (applied for denoising of sequencing error based on cluster generation by reads with only one nucleotide difference). Subsequently, taxonomic affiliation of DNA reads was performed using Rsyst::diatom barcoding library v4 (Rimet et al., 2016) and the naïve Bayesian method (Wang et al., 2007) with a confidence score threshold of 85%. All reads assigned to groups different from Bacillariophyta (diatoms) were excluded from further analyses. Similarity distance matrix based on pairwise distances between aligned reads (algorithm proposed by Needleman and Wunsch, 1970) was computed and used to cluster reads in operational taxonomic units (OTUs) with the furthest neighbor algorithm (95% similarity level as proposed by Mangot et al., 2013). Singleton removal and sample size normalization (to the smallest read abundance obtained among all samples) were performed and OTUs were taxonomically assigned based on consensus taxonomy of reads using confidence threshold of 80%. Taxonomy of abundant OTUs (above 100 reads), that were not precisely affiliated with the naïve Bayesian method, were verified with Rsyst::diatom blast to confirm their belonging to Bacillariophyta.

### Data Analysis

Quantitative (normalized abundances) and qualitative (presence–absence) matrices of OTUs were used to analyze diatom communities. All statistical analyses were performed using the free and open source R software (3.4.4, R development core team).

### Diatom Richness, Diversity, and Community Structure

Chao1 richness and Simpson diversity were calculated for all sampling dates and locations. NMDS was performed with Mothur Software on Bray–Curtis dissimilarity matrix to compare diatom communities. Differences between sites and seasons were tested with PERMANOVA (p < 0.05) (Anderson, 2001; McArdle and Anderson, 2001) made using the vegan package (Oksanen et al., 2017).

### Community Changes Linked to Environmental Factors

Redundancy analysis (RDA) was performed to infer the relationship between environmental factors and diatom communities from urban (U) and hospital (H) TE sites (vegan package, Oksanen et al., 2017). Prior to the analysis, the biological data (normalized quantitative OTU matrices) were Hellinger-transformed as recommended for linear ordination methods (Legendre and Gallagher, 2001). To reduce the impact of high correlations between environmental variables, an approach based on principal component analysis (PCA) was applied. Three pairs of highly correlated variables (r > 0.7) were identified: air temperature and irradiance (Temp.Irr), beta-blockers and NSAIDs (BB.NS), antibiotics and orthophosphate (ATB.PO4) and three separate PCAs were performed for each of the pairs. The first axis of each PCA (accounting for 98%, 92%, and 89% of the total variability for Temp.Irr, BB.NS, and ATB.PO4, respectively) was used in subsequent analysis as a synthetic variable representative for the variability of both variables.

In a next step, forward selection procedure was ran on the RDA performed on OTUs and the environmental variables ammonium, nitrite, nitrate, anticonvulsants, analgesics, TSS, COD, BB.NS, ATB.PO4, and Temp.Irr to obtain the most parsimonious RDA model with a reduced number of variables. The effect of time dependence was thereby considered by including the sampling month as a covariate. RDA was performed with the selected parameters, VIFs were verified and the significance of each variable and of the whole model was tested with permutation tests (999 permutations, p < 0.05).

Finally, variation partitioning based on multiple partial RDAs (pRDAs) was performed to quantify and test the proportion of biological variability explained by each of the selected variables (Borcard et al., 1992). The variation explained by each fraction was reported using unbiased adjusted R 2 (Peres-Neto et al., 2006). Permutation test (999 permutations, p < 0.05) was used to test the significance of each fraction.

### Relative Genera Abundances and Ecological Guild Classes

Relative abundances of OTUs affiliated on genus level and of ecological guild classes were calculated from means of sample triplicates and presented in barplots for each sampling location and period to study spatial and temporal diatom community dynamics. Ecological guilds were based on the original classification by Passy (2007) including the modifications proposed by Rimet and Bouchez (2012).

### Phylogeny and Indicator Species Analysis of OTUs

A phylogeny of sequences from the final list of OTUs was reconstructed. First, the OTU sequences were aligned against the rbcL sequences of the R-Syst::diatom database using the muscle algorithm (Edgar, 2004). Second, OTUs were inserted in a reference phylogeny of 604 species with the evolutionary placement algorithm (Berger et al., 2011) implemented in RAxML (Stamatakis, 2014). Finally, the tree was dated in relative time using PATHd8 (Britton et al., 2007).

Indicator species analysis was performed to identify OTUs that discriminate diatom communities in U, H, and RU locations. Indicator values based on relative abundance of OTUs were calculated following the method proposed by Dufrene and Legendre (1997) implemented in the labdsv package (Roberts, 2007). The indicator values of each OTU were calculated based on abundance values and relative frequency of occurrence and tested using a randomization procedure (999 permutations, p < 0.05). Each significant OTU was considered to be indicator of the location for which it exhibited the highest indicator value. Thus, a list of indicator OTUs was obtained for each location. Subsequently, indicator OTUs were denoted in the phylogenetic tree with location-specific colors for U (blue), H (black), and RU (green). Remaining non-significant OTUs were denoted with gray. Significant indicator OTUs of all three locations that appeared in RD were also labeled with the respective location-specific colors. Proportion of U, H, and RU indicator OTUs in RD communities was calculated (considering number of OTUs and DNA read abundances) to evaluate the influence of each location-characteristics on RD communities.

### Biological Diatom Index

Biological diatom index is a biotic index calculated after morphological identification of diatom species present in natural biofilms. It is routinely applied for water quality assessment of running waters in France and is used in the context of the WFD (Coste et al., 2009). In this study, OTUs were affiliated on species (or genus when species not available) level and their normalized quantitative matrices, based on read numbers, were used to calculate molecular BDI for each sampling location and period and to evaluate thereby the effect of TEs on water quality in the recipient river. BDI calculation was done with OMNIDIA 5 software (Lecointe et al., 1993).

### RESULTS

### Habitat Description

Principle major differences regarding nutrients and pharmaceuticals were found between sampling locations, especially distinguishing TE sites from river sites. Nutrients and PhC were clearly higher concentrated in the TEs (**Table 1**). Comparing urban and hospital effluents, the clearest trends were observed for orthophosphate and antibiotics with significantly higher concentrations in H, and beta-blockers, NSAIDs and ammonium – in U (paired Wilcoxon test, p < 0.05). In river locations, concentrations of PhC (especially antibiotics, NSAIDs and anticonvulsants) and orthophosphate showed rather increasing trend in RD compared to RU, which was not noticeable for ammonium, nitrite, nitrate, and TSS (Chonova et al., 2018a). During the observation period solar irradiance, air temperature and river flow increased from winter to summer. Threefold increase of the theoretical dilution factor for WWTP discharge was observed from winter to summer (**Table 2**).

### High-Throughput Sequencing

We obtained 3,760,984 reads from the HTS, with an average of 59,698 reads per sample and average length of 312 base pairs. After the bioinformatics filtration steps, 1,357,398 reads were retained and were clustered into 1,121 OTUs (**Supplementary**

TABLE 1 | Mean concentrations of nutrients and pharmaceuticals in urban (U) and hospital (H) treated effluents and Arve river up- (RU) and downstream (RD) from the WWTP Output (standard deviation in brackets).

PhC, pharmaceutical compounds; COD, chemical oxygen demand; TSS, total suspended solids; NSAIDs, non-steroidal anti-inflammatory drugs.

TABLE 2 | Mean values calculated from multiple daily measurements of solar irradiance, river temperature, precipitation, river flow, and WWTP discharge for each colonization period (standard deviation in brackets; data sources: INRA meteorological station<sup>∗</sup> , federal office of environment of Switzerland∗∗ and SIPIBEL observatory).

**Table S1**). All samples were subsampled to 11,306 reads (lowest read abundance obtained for a sample) to allow intersample comparison which resulted in total of 1,076 OTUs with minimum of 103 and maximum of 301 OTUs per sample. R-Syst::diatom v4 database enabled taxonomic affiliation of 103 different species from 48 different genera. 81% of the OTUs (90% of the reads) were affiliated at family level, 75% of the OTUs (85% of the reads) – at genus level and 57% of the OTUs (70% of the reads) – at species level. Those results are in consistence with other studies (e.g., Vasselon et al., 2017b). Most of the unclassified reads were found in H – a location that represented the most extreme environment in terms of concentrations of micro- and macropollutants. However, results from R-Syst::diatom blast confirmed that unclassified OTUs belong to the group Bacillariophyta and were not misclassified.

### Diatom Community Changes and the Role of Micro- and Macropollutants

Boxplots on **Figure 2** compare richness (Chao1) and diversity (Shannon) of benthic diatom communities between sampling locations and suggested increased richness and diversity in the river compared to the TEs. Higher richness and diversity were observed in H than in U. Comparing river locations, alpha diversity changed depending on the method – Chao1 reported higher richness in RU, and Shannon suggested tendency for higher diversity in RD. Seasonal trends of diatom richness and diversity were not observed.

Non-metric multidimensional scaling analysis of quantitative OTU matrix showed that diatom communities were mainly grouped by sampling location (**Figure 3A**). Significant differences between U and H communities were represented on the first axis, and between TE sites and river sites on the second axis (PERMANOVA, p < 0.05). RD communities were more similar to U and H than RU communities. As secondary shaping factor, seasonal changes were visible for all sampling locations.

Redundancy analysis was performed to explore the relationship between diatom community structure and environmental variables measured in TEs. Forward selection procedure identified that the best parsimonious model includes the variables ATB.PO4 (representing 89% of the common variability of antibiotics and orthophosphate) and BB.NS (representing 92% of the common variability of beta-blockers and NSAIDs) as the most important for diatom community changes (**Figure 3B**). RDA model including ATB.PO4 and BB.NS was significant and the two variables tested separately were significant to the model (p < 0.05). The low VIF values (VIF < 3) suggested that there is no issue with collinearity. The first two axes of the model accounted for 61% and 3% of the variability in diatom community structure, respectively, when considering the effect of time dependence as covariate. The first axis defined a local gradient between the two basins, where diatom community development in U was rather linked to higher concentrations of beta-blockers and NSAIDs and in H – to antibiotics and orthophosphate. The second axis explained a low proportion of the variance only. Variation partitioning revealed that the effect explained by ATB.PO4 and BB.NS separately remained significant (p < 0.05) and it accounted for 17% and 11% of the total variability for ATB.PO4 and BB.NS, respectively. The variance shared between the two variables was 36%.

### Spatial and Temporal Variations in Diatom Genera and Ecological Guild Classes

Relative abundances of diatom genera (a) and ecological guild classes (b) per sampling site and period are presented in **Figure 4**. Planktic species were rare in the biofilms. Clear dichotomy in ecological guild classes was found between TE sites and river sites.

At TE sites mainly the development of polysaprobic species (e.g., Sellaphora, Craticula) was observed. Motile groups (represented mainly by the genera Craticula, Sellaphora, and Nitzschia) dominated in treated wastewater with relative abundance up to 95% in U. Their dominance slightly decreased in H, where motile guild was partially replaced by low- (Achnanthidium and Planothidium) and high-profile groups (Gomphonema). Seasonal trends were observed in U for the genus Craticula and Sellaphora, as Craticula was dominating in colder months (February and March) with ca. 70% relative abundance but was replaced by Sellaphora in warmer months. In H, Nitzschia, Gomphonema, and Fragilaria exhibited slightly higher relative abundance in colder months and decreased in summer.

In river communities, genera relative abundances were more evenly distributed and low- and high-profile species were better

represented (especially in RU). Importance of Craticula and Sellaphora decreased and motile groups were rather represented by Mayamaea, Navicula, and Nitzschia. Overall, in the river, groups with oligotrophic and oligosaprobic (e.g., Achnanthidium, Encyonema) preferences were favored. This trend was clearly stronger expressed in the river upstream from the WWTP output. Seasonal trends were observed for low- and high-profile diatoms. Low-profile groups (mainly Achnanthidium, Cocconeis, and Reimeria) were more abundant in warmer months, while increase in high-profile diatoms (mainly Encyonema and Gomphonema) was observed rather in colder months, especially upstream from the WWTP.

### Location-Specific Indicator OTUs

Abundances of OTUs were calculated as mean of all samples for each location and are presented (log-transformed) in **Figure 5**. Indicator OTUs of U, H, and RU (represented in blue, black, and green on the figure, respectively) were located in

different parts of the phylogenetic tree. Indicator OTUs of U belonged mainly to Nitzschia, Sellaphora, and Craticula, whereas a large phylogenetic clade comprising Cocconeis, Achnanthidium, Planothidium, Encyonema, and Gomphonema was hardly present in this location. By contrast, many indicator OTUs belonging mainly to Achnanthidium, Planothidium, and Gomphonema appeared again in H, and became even more important for RU. Indicator OTUs of RU were detected throughout the phylogenetic tree, except for Craticula and Sellaphora. These results confirmed general trends revealed with the relative abundance of genera.

**Figure 5** clearly illustrates that most of the OTUs found in RD corresponded to indicator OTUs of RU (53% of the OTUs, corresponding to 88% of the DNA reads). However, numerous indicator OTUs of U or H that were not present in RU appeared in RD (e.g., Craticula, Sellaphora, and Nitzschia). This proportion was 14% of the OTUs (corresponding to 7% of the reads) for U, and 13% of the OTUs (corresponding to 2% of the reads) for H.

# Biological Diatom Water Quality Index

Biological diatom index was calculated from molecular data and represented in **Figure 6**. As expected, BDI suggested that water quality was lower at TE sites compared to the recipient river. Lowest values were observed for H (quality between "bad" and "poor"), followed by U (between "poor" and "moderate"). At river sites, BDI revealed "high" quality status in RU and decrease from "high" to "good" in RD.

### DISCUSSION

The development and behavior of diatom communities in response to TEs remain poorly documented. Microalgae have been previously associated to wastewater treatment due to their ability to assimilate nitrogen and phosphorus and to bind heavy metals (e.g., Congestri et al., 2005). In this

context, in recent studies Ghosh and Love (2011) investigated algal assemblages from the water column of a secondary WWTP tank, and Congestri et al. (2005) reported communities developed on tank walls and artificial substrata in WWTP sedimentation tanks. In such highly polluted environments, diatoms represent one of the most abundant groups among algae (e.g., Ghosh and Love, 2011). Although diatoms are regularly used for biomonitoring and their sensitivity to micropollutants is well known, species composition and structural dynamics in WWTPs and their TEs have rarely been documented in this context. The present study reports for the first time the reliability of DNA metabarcoding to detect spatial and temporal

responses of diatom assemblages to high environmental pressure in urban and hospital TEs and contributes to our understanding of the effluents' impact on community dynamics in the recipient aquatic environment.

Parallel analysis of nutrients, PhC and seasonal parameters provided essential information about the characterization of different habitats and confirmed higher presence of pollutants in the TEs compared to the river. As expected, comparison of U and H showed that hospital TEs contained higher total concentrations of PhC (especially antibiotics) and phosphate caused by the regular use of PhC and specific detergents in the hospital. However, better removal efficiency of NSAIDs and beta-blockers in the hospital treatment process might have resulted in their higher concentrations in urban TEs (Chonova et al., 2016). Release of WWTP effluents in the recipient river led to increased concentrations of PhC (especially NSAIDs, antibiotics, and anticonvulsants) and phosphate downstream from the TEs output. In contrast nutrients such as ammonium, nitrate, nitrite, and TSS did not show such up – to downstream gradient. Temperature, solar irradiance and river flow (caused by the melting of glacier) followed the usual seasonal gradient with increase in summer. No extreme water level changes or flood events were observed during the study experimental period.

# Composition of Diatom Communities Linked to Habitats Characteristics

Diversity and community structure of benthic diatoms exhibited clearly stronger local changes than seasonal ones, which highlights the relevance of location-specific factors (e.g., nutrients and PhC) over seasonal factors (e.g., temperature and solar irradiance) (**Figures 2**, **3**). The high relevance of nutrients in structuring diatom communities has been well studied (e.g., Patrick, 1961; Lange-Bertalot, 1979; Larson and Passy, 2012; Marcel et al., 2013). Several studies also show significant effect of PhC (e.g., NSAIDs, beta-blockers, antibiotics) and other micropollutants (pesticides, heavy metals, etc.) on benthic diversity and community structure (e.g., Hagenbuch and Pinckney, 2012; Larras et al., 2013; Corcoll et al., 2014). Temperature, solar irradiance and flow velocity (for river sites) were previously shown to be crucial factors for the development of benthic diatom communities (Lange et al., 2011; Villeneuve et al., 2011; Larras et al., 2013). However, in this study they remained of secondary importance. Such primary local and secondary seasonal community dynamics at the same study site were reported for bacterial biofilms after short- and long-term biofilm colonization periods, where both nutrients and PhC played a key role in the shaping of the community structure (Chonova et al., 2016, 2018a).

### Response of Diatoms to Pollutants in Urban and Hospital Treated Effluents

Communities developed at TE sites and at river sites differed strongly regarding their richness and diversity which reflected the contrasting ecological differences between the two systems (**Table 1**). Considerably lower diatom richness and diversity were observed at TE sites which were the most severely

disturbed and showed higher concentrations of micro- and macropollutants compared to the natural river system (**Figure 2**). Such low diatom diversity was also reported in studies investigating algal communities in WWTP clarification tanks and is probably linked to the high anthropogenic pressure which leads to the elimination of less resistant species (e.g., Sládecková et al., 1983 ˇ ; Congestri et al., 2005). However, restricted colonization possibilities in these "artificial systems" offering limited microhabitats compared to natural environments may also affect the number of species. This question needs to be further investigated applying translocation experiments as described for example by Proia et al. (2013b).

Diatom community structure at TE sites also differed from river sites regarding HTS-OTUs, genera taxonomic level and ecological guild classes (**Figures 3**, **4**). These trends allowed to better understand the community dynamics linked to the ecological preferences of diatoms, and clearly showed that community adaptations are in accordance with habitat characteristics. In general, TE sites characterized by higher anthropogenic pressure and lower water velocity and turbulence, were dominated by few genera with polysaprobic preferences and well-developed raphe system belonging to the motile guild (e.g., Craticula, Nitzschia, Sellaphora, etc.). Some of these genera (Craticula, Nitzschia, Sellaphora, and Gomphonema) were also previously reported as dominant in WWTP tanks (Congestri et al., 2005). As described by Passy (2007) and Rimet and Bouchez (2012), the motile guild included free-moving species with tolerance to high nutrient concentrations and low resistance to flow velocity. Higher tolerance of motile species to micropollutants (e.g., herbicides and fungicides) has also been reported previously (Rimet and Bouchez, 2011). This resistance to pollution is linked to the ability of motile diatoms to optimize their position in the biofilm and thereby avoid disturbances (Lengyel et al., 2015). The adaptations of motile groups to survive and reproduce in highly contaminated environments explains their dominance over other guilds in TEs. Dominance of raphid taxa in the TEs was also reported previously in WWTPs (Congestri et al., 2005; Ghosh and Love, 2011) and may be explained by their higher resistance to micropollutants which may lead to the replacement of higher sensitive centric and araphid taxa (e.g., Larras et al., 2014 for pesticides).

According to Ghosh and Love (2011), algal diversity and composition may differ strongly between different WWTPs. Here, we also found differences when comparing separately treated urban and hospital effluents, despite their close geographic proximity. Hence, these differences are probably rather linked to the origin-specific effluent composition than to the spatial effect on colonization possibilities. Interestingly, lower richness and diversity was observed in U, where concentrations of total nutrients and PhC were lower (**Figure 2** and **Table 1**). Differences between U and H communities were also observed on OTU and genera level, and to lower extent – regarding functional traits (**Figures 3**, **4**). Urban communities in colder months (February and March) were dominated by Craticula, which was replaced by Sellaphora in warmer months. Both species are symmetrical biraphid diatoms and exhibit similar trophic preferences since they tolerate nutrient rich environments. Their similar characteristics may lead to competition resulting in suppression of Craticula in warmer – and Sellaphora in colder months. Nitzschia was also regularly found in U, but did not exhibit such strong seasonal variations. In hospital communities, low- (Achnanthidium and Planothidium) and high-profile groups (Gomphonema) showed higher relative abundances compared to U (**Figure 4B**). However, motile groups (Craticula, Nitzschia, and Sellaphora) remained dominant. Regarding seasonal trends, Nitzschia, Gomphonema, and Fragilaria exhibited slightly higher relative abundance in colder months, but genera composition in H remained generally stable between seasons.

Aiming to better comprehend community changes and to categorize the importance of nutrients, pharmaceuticals, and seasonal conditions, diatom community dynamics in U and H were further studied on OTU level with RDA. Forward selection procedure defined beta-blockers, NSAIDs, antibiotics and phosphate as the most important factors explaining a significant part of the diatom OTU variability (64%). Higher concentrations of beta-blockers and NSAIDs (accounting for 11% of the biological variability) were rather linked to urban communities, while antibiotics and phosphate (accounting for 17%) – to hospital ones. Beta-blockers and NSAIDs, in the same concentration range as reported in U, may inhibit algal photosynthesis processes (Bonnineau et al., 2010) and change benthic community structure (Corcoll et al., 2014). Lower diversity in U may therefore be linked to the influence of these PhC. Other studies showed that NSAIDs may change phosphatase activity (Proia et al., 2013b) and reduce photosynthesis (Ding et al., 2017). Negative effect of antibiotics on benthic communities has also been reported (e.g., Proia et al., 2013b) and associated with reduction or elimination of diatom mobility (Hagenbuch and Pinckney, 2012) and inhibition of diatoms growth due to reduced photosynthesis (Pinckney et al., 2013; Guo et al., 2016). Hence, mobility reduction in H, where antibiotic concentrations were more than 20 times higher than in U, may have led to discrimination of motile diatoms, favoring development of low- and high-profile species. Furthermore, effect of antibiotics on bacterial communities may lead to indirect changes in diatom assemblages (e.g., Windler et al., 2015). The remaining PhC (paracetamol and carbamazepine) were excluded by the forward selection procedure in the course of RDA. However, considering the high concentrations of the analgesic paracetamol and the anticonvulsant carbamazepine, we cannot exclude their possible effect on diatom communities. Indeed, paracetamol may also impact photosynthesis and cause algal growth inhibition (Proia et al., 2013b). The presence of non-measured PhC, metabolites, pesticides, heavy metals and detergents may also have played a role in the shaping of diatom communities (e.g., Morin et al., 2009; Ricciardi et al., 2009). The mixture of highly concentrated pollutants found in threated effluents may exhibit additive, synergistic and antagonistic effects on benthic communities (Hagenbuch and Pinckney, 2012; Verlicchi, 2018). Such complex cocktail effects may strongly vary depending on the present compounds and the microbial communities. Hence, single-compound/species studies are not sufficient to make meaningful predictions

of the ecological impact of PhC and potential toxic effects of "micropollutant cocktails" need to be better understood (Hagenbuch and Pinckney, 2012).

### Dynamics of Benthic Diatoms in Recipient River

River habitats (RU and RD) in this study are characterized as natural aquatic ecosystem and exhibit lower anthropogenic pressure and higher flow velocity and turbidity (particularly in spring and summer when river flow increases due to the snow and glacier melting) compared to TE sites. In the river, species with oligosaprobic preferences were generally favored and the development of diverse diatom communities and higher evenness between genera was facilitated. Consequently, lower dominance of motile groups (here represented mainly by Mayamaea, Navicula, and Nitzschia) and better growth of low- and high-profile groups was observed. In both river locations, higher prevalence of low-profile genera (represented mainly by Achnanthidium, Cocconeis, and Reimeria) was observed in summer (June and July). This can be explained by their better adaptation to resist to high current velocity and strong water turbulence that increase in these months. In contrast, high-profile groups (represented mainly by Encyonema and Gomphonema), that are not well adapted to high velocity, developed better in winter (Passy, 2007).

Benthic diatom assemblages in the river up- and downstream from the WWTP effluent were compared to evaluate the effect of TEs on the recipient aquatic environment. Regarding richness and diversity, we observed inconsistency between Chao1 (lower richness in RD) and Shannon (lower diversity in RU) (**Figure 2**). This difference may be linked to the increase in anthropogenic pressure in RD resulting in loss of rare species (ca. 20% less OTUs with <3 reads) (Proia et al., 2013a). Biggs (2000) observed that under certain circumstances, major species replacements with increasing eutrophication may also result in increase of species richness. Contrasting conclusions have been reported for the response of diatom richness and diversity to presence of micropollutants as rather a decrease in the presence of micropollutants (mainly herbicides and heavy metals) was observed (e.g., Genter and Lehman, 2000; Sabater, 2000; Morin et al., 2009; Ricciardi et al., 2009), while others were unable to detect such relation (e.g., Hirst et al., 2002; Marcel et al., 2013). Richness and diversity indices are not always relevant to assess waterbodies pollution levels (Ricciardi et al., 2009; Pandey et al., 2017). They reduce community information to a single number, which leads to huge loss of essential information about various crucial community characteristics (e.g., taxonomic structure, functional traits, phylogenetic relations, etc.) (Biggs, 2000).