Teofana Chonova1,2*
Rainer Kurmayer1
Frédéric Rimet2Jérôme Labanowski3Valentin Vasselon2François Keck2,4
Paul Illmer5
Agnès Bouchez2- 1Research Department for Limnology, Mondsee, Faculty of Biology, University of Innsbruck, Mondsee, Austria
- 2UMR CARRTEL, INRA, Université Savoie Mont Blanc, Thonon-les-Bains, France
- 3UMR IC2MP 7285, CNRS, Université de Poitiers, ENSIP, Poitiers, France
- 4Department of Aquatic Sciences and Assessment, Faculty of Natural Resources and Agricultural Sciences, Swedish University of Agricultural Sciences, Uppsala, Sweden
- 5Department of Microbiology, Faculty of Biology, University of Innsbruck, Innsbruck, Austria
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.
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, cost-efficient 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).
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.
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 cm2).
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 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 cm2) 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 cm2 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 (m3.d) by the average WWTP discharge (m3.d−1).
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® polymerase and the primer pairDiat_rbcL_708F (AGGTGAAGTTAAAGGTTCATACTTDAA) (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® 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® 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 R2 (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).
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).
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∗, federal office of environment of Switzerland∗∗ and SIPIBEL observatory).
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 S1). All samples were subsampled to 11,306 reads (lowest read abundance obtained for a sample) to allow inter-sample 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.
Figure 2. Boxplots presenting variation in (A) Chao1 richness and (B) Shannon diversity calculated from HTS-OTUs of benthic diatoms developed in urban (U) and hospital (H) treated effluent sites and river sites up- (RU) and downstream (RD) from the WWTP output.
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.
Figure 3. (A) NMDS two-dimensional plot of observed similarities between rbcL OTU profiles of benthic diatom communities of urban (U) and hospital (H) treated effluents and river sites up- (RU) and downstream (RD) from the WWTP output (stress value = 0.22). (B) Biplot from redundancy analysis based on diatom OTU profiles for periphytic samples from TE sites (U and H), and the variables ATB.PO4 and BB.NS chosen with forward selection procedure. Colors delineate sampling locations: U (in blue), H (in black or gray), RU (in green) and RD (in red). Numbers (1, 2, and 3) in the sample name denote triplicates.
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.
Figure 4. Relative abundances of diatom (A) genera and (B) ecological guild classes (in %) in each sampling site and period. “Others” – sum of all genera representing less than 5%.
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. Phylogenetic tree of 1,076 OTUs (log-transformed relative average abundance) of benthic diatoms from urban (U) and hospital (H) treated effluent sites and river sites upstream (RU) and downstream (RD) from the WWTP output. Colors delineate indicator OTUs for each location: U (in blue), H (in black), and RU (in green).
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.
Figure 6. Biological diatom water quality index (BDI) calculated for benthic diatoms from urban (U) and hospital (H) treated effluent sites and river sites up- (RU) and downstream (RD) from the WWTP.
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ádečková 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).
Changes in community structure between river sites reflected the pollution induced by the TEs discharge that resulted in increasing trends of phosphate, NSAIDs, antibiotics and anticonvulsants downstream from the WWTP output. Species with oligotrophic and oligosaprobic preferences and high-profile taxa (that are more sensitive to contamination of micropollutants, Marcel et al., 2013) were rather disadvantaged downstream from the WWTP comparing to upstream. In RD, these were replaced by motile groups, reflecting the increase in anthropogenic pollution (Rimet and Bouchez, 2011). Motile groups are usually less resistant to current velocity (Passy, 2007) and it was previously reported that they develop better in winter (Stenger-Kovács et al., 2013). However, in this study they exhibited such seasonal trend in RD only, which may be rather linked to the lower dilution factor of TEs in winter, potentially leading to stronger impact on river diatoms during winter.
Chronic input of trace emerging contaminants may in general lead to decrease in biofilm biomass and primary productivity and may thereby influence the ecology of benthic diatom assemblages (Hagenbuch and Pinckney, 2012). Such changes in more sensitive systems may impair higher trophic levels and lead to important alterations in river ecosystem functioning (Pinckney et al., 2013). Long-term exposure to contaminants may finally result in more tolerant assemblages (e.g., Corcoll et al., 2014). Results presented here are helpful to follow the changes on natural diatom communities potentially caused by the release of WWTP effluents. However, further studies are needed to better comprehend the cocktail effect of PhC and other micropollutants on benthic diatoms in aquatic environments. When studying such relationships, the potential effect of confounding factors (e.g., eutrophication levels, seasonal conditions) has to be carefully considered (Marcel et al., 2013).
Tracking the Effect of Urban and Hospital Treated Effluents on River Communities
Indicator Species Analysis
Indicator species analysis in combination with phylogenetic tree representation was helpful to better understand the effect of TEs on diatom assemblages and disentangle community changes in natural aquatic environment receiving urban and hospital TEs. Clades of phylogenetically close indicator OTUs behaved similarly and appeared simultaneously in the different habitats. Such patterns were expected from highly similar OTUs, as the clustering method applied here may generate several OTUs belonging to the same species. However, this trend also appeared on a larger scale in the phylogeny. Indicator OTUs from a large clade including mainly Cocconeis, Achnanthidium, Planothidium, Encyonema, and Gomphonema were hardly present in U, and those from a clade including mainly Navicula and Mayamaea were not found in H. Such phylogenetically related groups are likely to share similar characteristics like guild classes and ecological preferences (Keck et al., 2016). Cocconeis, Achnanthidium, Planothidium, Encyonema, and Gomphonema all possess a developed raphid system and belonged to the high- or low-profile guilds. In contrast, Navicula and Mayamaea are part of the motile guild. Indicator OTUs discriminating RU were detected throughout the phylogenetic tree. Entire clades from genera highly dominant at TE sites (e.g., Craticula, Sellaphora, and partially Nitzschia) were hardly found in RU implying their ecological optimum in contaminated environments.
According to the ISA, RD communities showed highest similarity to RU communities and were composed by 53% of OTUs (88% of DNA reads) indicative of RU. However, the two river locations differed significantly. This can be explained by the continuous influence of micro- and macropollutants released by TEs that may inhibit the growth of sensitive species (Proia et al., 2013a), but also by possible species colonization downstream. Furthermore, 27% of the RD OTUs (9% of DNA reads) corresponded to indicator OTUs of U or H communities. Part of these OTUs probably corresponds to species that were released with the TEs and captured by the natural river biofilm downstream from the WWTP output. Motile, high-profile and especially planktic groups are more likely to disperse downstream with the water flow due to their lower resistance to current velocity (Liu et al., 2013; Dong et al., 2016). Indicator OTUs discriminating U or H and belonging to Fragilaria, Gomphonema (high-profile), and Cyclotella (planktic) were found in RD communities, despite their absence in RU. Transportation of free DNA from dead cells by the water flow may also occur and this DNA may be captured in river biofilms (e.g., Deiner et al., 2016; Pont et al., 2018). More directly, the permanent release of TEs causing constant pressure (PhC and nutrient load) in RD may favor the maintenance of more tolerant species and promote their stabilization in the biofilm.
River downstream communities seemed to be more influenced by U than by H communities, as 7% of the OTU reads found in RD corresponded to indicator OTUs of U, and only 2% – to indicator OTUs of H. The urban WWTP discharge was more than 15 times higher than the hospital one, which probably led to higher influence on water chemistry changes in RD.
Biological Diatom Index
Finally, molecular BDI was calculated from taxonomically assigned OTUs to evaluate the effect of TEs’ release on water quality in the recipient river. The application of BDI is usually limited to natural aquatic environments. However, WFD stations may be located in channels in urban environment build up with substrate different from natural river (similar to our case), and pollutant concentrations measured at TE sites during this study are comparable to these reported for highly polluted natural systems. Hence, we calculated BDI for all habitats (including TE sites) to evaluate the response of the entire diatom community to high anthropogenic pressure and track the effect of TEs on river communities. Compared to the river sites, communities in U and H reflected relatively low water quality (between “bad” and “moderate”), responding to the high anthropogenic pressure. WWTP discharges are generally monitored to avoid large alterations in river quality. The studied WWTP respected all national and European norms linked to removal efficiency and pollutants release (Chonova et al., 2016). Nevertheless, BDI based on molecular data suggested degradation of the water quality downstream from the WWTP output. This finding was confirmed by classical BDI based on microscopy counts (Sipibel Report, 2016) and it pointed out the efficiency of BDI to highlight changes in river ecology linked to the release of TEs.
Conclusion
The present study shows that assessment of benthic diatoms using DNA metabarcoding is efficient to detect spatial and temporal community responses to pharmaceutical pressure in TEs of urban and hospital wastewaters and contributes to our understanding of the potential effluent impact on community dynamics when released in the recipient aquatic environment. We detected habitat-specific changes linked to the WWTP effluents on different community levels – richness, diversity, taxonomic composition, functional traits, and phylogenetic position. The changes in taxonomy and ecological traits included shift in proportion of polysaprobic motile groups in TE vs. oligosaprobic/oligotrophic groups belonging to the low- and high-profile guild in the river (especially upstream). RDA suggested that beta-blockers, NSAIDs, antibiotics and phosphate were among the most important factors driving community dynamics in TEs. ISA was helpful to evaluate the negative effect of these TEs on natural river community composition, revealing that 27% of OTUs detected in RD communities were indicative of urban or hospital treated effluent origin. Those OTUs may be either directly transferred with the TE and captured in the biofilm or maintained and developed further in RD due to the continuous release of TEs that changes chemical water parameters. Finally, the BDI calculated to evaluate the ecological status of the recipient river suggested that the release of TEs may lead to water quality decrease.
From these results, we can conclude that there is clearly a change in community structure linked to the WWTP effluents. As discussed above, studying environmental relationships between microbial communities and pollutants remains challenging due to occurring intercorrelations and the potential effect of confounding factors. Study on larger scale in combination with laboratory experiments would be helpful to explore and strengthen further relations between water chemistry and diatom community changes and to confirm findings reported here.
Author Contributions
JL, AB, and TC contributed to the experimental design and sampling. TC contributed to the molecular experiments and data collection and organization. VV and TC contributed to the bioinformatics. RK, FK, FR, PI, and TC contributed to data analysis and statistics. RK, FR, JL, AB, and TC contributed to the interpretation of the data. FK and TC contributed to the preparation of the figures. RK, FR, JL, VV, FK, PI, AB, and TC contributed to the revision and approval of the manuscript.
Funding
This study was done as part of the SIPIBEL field observatory on hospital effluents and urban WWTPs and was funded by ANSES (French National Agency for Food, Environmental and Occupational Health and Safety), project “PERSIST-ENV” (Effluents hospitaliers et persistance environnementale des médicaments et de bactéries pathogènes) #2012/2/149 as part of the programme Environnement-Santé-Travail (French Ministers in charge of ecological and environmental issues). An Excellence Scholarship for doctoral programmes was obtained from the University of Innsbruck, Office of the Vice Rector for Research.
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.
Acknowledgments
The article is based on collaborative work with COST Action DNAqua-Net (CA15219), supported by the COST (European Cooperation in Science and Technology) program. We thank Frankreich-Schwerpunkt (University of Innsbruck), and the EU co-funded project Interreg Alpine Space “Eco-AlpsWater” for their support. We thank L. Wiest (ISA, Lyon) and V. Lecomte (GRAIE) for their contribution.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2019.00653/full#supplementary-material
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, non-steroidal 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.
References
AFNOR (1997). French Standard Operating Procedures, Water Quality-Analytical Methods, 2nd Edn, Vol. 3. La Défense: AFNOR.
Anderson, M. J. (2001). A new method for non-parametric multivariate analysis of variance. Austral Ecol. 26, 32–46.
Berger, S. A., Krompass, D., and Stamatakis, A. (2011). Performance, accuracy, and web server for evolutionary placement of short sequence reads under maximum likelihood. Syst. Biol. 60, 291–302. doi: 10.1093/sysbio/syr010
Biggs, B. J. F. (2000). New Zealand Periphyton Guideline: Detecting, Monitoring and Managing Enrichment of Streams. Wellington: Ministry for the Environment, 121.
Bonnineau, C., Guasch, H., Proia, L., Ricart, M., Geiszinger, A., Romaní, A. M., et al. (2010). Fluvial biofilms: a pertinent tool to assess b-blockers toxicity. Aquat. Toxicol. 96, 225–233. doi: 10.1016/j.aquatox.2009.10.024
Borcard, D., Legendre, P., and Drapeau, P. (1992). Partialling out the spatial component of ecological variation. Ecology 73, 1045–1055. doi: 10.2307/1940179
Britton, T., Anderson, C. L., Jacquet, D., Lundqvist, S., and Bremer, K. (2007). Estimating divergence times in large phylogenetic trees. Syst. Biol. 56, 741–752. doi: 10.1080/10635150701613783
Bruder, K., and Medlin, L. K. (2007). Molecular assessment of phylogenetic relationships in selected species/genera in the naviculoid diatoms (Bacillariophyta). I. The genus Placoneis. Nova Hedwig. 85, 331–352. doi: 10.1127/0029-5035/2007/0085-0331
Chonova, T., Keck, F., Labanowski, J., Montuelle, B., Rimet, F., and Bouchez, A. (2016). Separate treatment of hospital and urban wastewaters: a real scale comparison of effluents and their effect on microbial communities. Sci. Total Environ. 542, 965–975. doi: 10.1016/j.scitotenv.2015.10.161
Chonova, T., Labanowski, J., and Bouchez, A. (2017). “Contribution of hospital effluents to the load of micropollutants in WWTP influents,” in Hospital Wastewaters – Characteristics, Management, Treatment and Environmental Risks. The Handbook of Environmental Chemistry, ed. P. Verlicchi (Heidelberg: Springer), 135–152. doi: 10.1007/698_2017_21
Chonova, T., Labanowski, J., Cournoyer, B., Chardon, C., Keck, F., Laurent, E., et al. (2018a). River biofilm community changes related to pharmaceutical loads emitted by a wastewater treatment plant. Environ. Sci. Pollut. Res. Int. 25, 9254–9264. doi: 10.1007/s11356-017-0024-0
Chonova, T., Lecomte, V., Krajewski, J. L., Bouchez, A., Labanowski, J., Dagot, C., et al. (2018b). The SIPIBEL project: treatment of hospital and urban wastewater in a conventional urban wastewater treatment plant. Environ. Sci. Pollut. Res. Int. 25, 9197–9206. doi: 10.1007/s11356-017-9302-0
Congestri, R., Cox, E. J., Cavacini, P., and Albertano, P. (2005). Diatoms (Bacillariophyta) in phototrophic biofilms colonising an Italian wastewater treatment plant. Diatom Res. 20, 241–255. doi: 10.1080/0269249X.2005.9705634
Corcoll, N., Acuña, V., Barceló, D., Casellas, M., Guasch, H., Huerta, B., et al. (2014). Pollution-induced community tolerance to non-steroidal anti-inflammatory drugs (NSAIDs) in fluvial biofilm communities affected by WWTP effluents. Chemosphere 112, 185–193. doi: 10.1016/j.chemosphere.2014.03.128
Coste, M., Boutry, S., Tison-Rosebery, J., and Delmas, F. (2009). Improvements of the Biological Diatom Index (BDI): description and efficiency of the new version (BDI-2006). Ecol. Indic. 9, 621–650. doi: 10.1016/j.ecolind.2008.06.003
Deiner, K., Fronhofer, E. A., Mächler, E., Walser, J.-C., and Altermatt, F. (2016). Environmental DNA reveals that rivers are conveyer belts of biodiversity information. Nat. Commun. 7:12544. doi: 10.1038/ncomms12544
Ding, T., Yang, M., Zhang, J., Yang, B., Lin, K., Li, J., et al. (2017). Toxicity, degradation and metabolic fate of ibuprofen on freshwater diatom Navicula sp. J. Hazard. Mater. 330, 127–134. doi: 10.1016/j.jhazmat.2017.02.004
Dong, X., Li, B., He, F., Gu, Y., Sun, M., Zhang, H., et al. (2016). Flow directionality, mountain barriers and functional traits determine diatom metacommunity structuring of high mountain streams. Sci. Rep. 6:24711. doi: 10.1038/srep24711
Dufrene, M., and Legendre, P. (1997). Species assemblages and indicator species: the need for a flexible asymmetrical approach. Ecol. Monogr. 67, 345–366. doi: 10.2307/2963459
Edgar, R. C. (2004). MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792–1797. doi: 10.1093/nar/gkh340
Edgar, R. C., Haas, B. J., Clemente, J. C., Quince, C., and Knight, R. (2011). UCHIME improves sensitivity and speed of chimera detection. Bioinformatics 27, 2194–2200. doi: 10.1093/bioinformatics/btr381
Esteves, S. M., Keck, F., Almeida, S. F. P., Figueira, E., Bouchez, A., and Rimet, F. (2017). Can we predict diatoms herbicide sensitivities with phylogeny? Influence of intraspecific and interspecific variability. Ecotoxicology 26,1065–1077. doi: 10.1007/s10646-017-1834-z
European Commission (2013). Directive 2013/39/EU of the European Parliament and of the Council of 12 August 2013 Amending Directives 2000/60/EC and 2008/105/EC as Regards Priority Substances in the Field of Water Policy (2011/0429 (COD)). Brussels: European Commission.
Genter, R. B., and Lehman, R. M. (2000). Metal toxicity inferred from algal population density, heterotrophic substrate use, and fatty acid profile in a small stream. Environ. Toxicol. Chem. 19, 869–878. doi: 10.1002/etc.5620190413
Ghosh, S., and Love, N. G. (2011). Application of RbcL based molecular diversity analysis to algae in wastewater treatment plants. Bioresour. Technol. 102,3619–3622. doi: 10.1016/j.biortech.2010.10.125
Guo, J., Selby, K., and Boxall, A. B. (2016). Assessment of the risks of mixtures of major use veterinary antibiotics in European surface waters. Environ. Sci. Technol. 50, 8282–8289. doi: 10.1021/acs.est.6b01649
Hagenbuch, I. M., and Pinckney, J. L. (2012). Toxic effect of the combined antibiotics ciprofloxacin, lincomycin, and tylosin on two species of marine diatoms. Water Res. 46, 5028–5036. doi: 10.1016/j.watres.2012.06.040
Hirst, H., Jüttner, I., and Ormerod, S. J. (2002). Comparing the responses of diatoms and macroinvertebrates to metals in upland streams of Wales and Cornwall. Freshw. Biol. 47, 1752–1765. doi: 10.1046/j.1365-2427.2002.00904.x
Keck, F., Bouchez, A., Franc, A., and Rimet, F. (2016). Linking phylogenetic similarity and pollution sensitivity to develop ecological assessment methods: a test with river diatoms. J. Appl. Ecol. 53, 856–864. doi: 10.1111/1365-2664.12624
Kermarrec, L., Franc, A., Rimet, F., Chaumeil, P., Humbert, F., and Bouchez, A. (2013). Next-generation sequencing to inventory taxonomic diversity in eukaryotic communities: a test for freshwater diatoms. Mol. Ecol. Resour. 13, 607–619. doi: 10.1111/1755-0998.12105
Labanowski, J., Laurent, E., Chonova, T., Bouchez, A., Cournoyer, B., Marjolet, L., et al. (2016). Hospital effluents and environ- mental persistence of pathogenic bacteria and pharmaceutical com- pounds – the Persist-Env approach. Tech. Sci. Méthodes 6, 22–30. doi: 10.1051/tsm/201606022
Lange, K., Liess, A., Piggott, J. J., Townsend, C. R., and Matthaei, C. D. (2011). Light, nutrients and grazing interact to determine stream diatom community composition and functional group structure. Freshw. Biol. 56, 264–278. doi: 10.1111/j.1365-2427.2010.02492.x
Lange-Bertalot, H. (1979). Pollution tolerance of diatoms as a criterion for water quality estimation. Nova Hedwig. 64, 285–304.
Larras, F., Keck, F., Montuelle, B., Rimet, F., and Bouchez, A. (2014). Linking diatom sensitivity to herbicides to phylogeny: a step forward for biomonitoring? Environ. Sci. Technol. 48, 1921–1930. doi: 10.1021/es4045105
Larras, F., Lambert, A.-S., Pesce, S., Rimet, F., Bouchez, A., and Montuelle, B. (2013). The effect of temperature and a herbicide mixture on freshwater periphytic algae. Ecotoxicol. Environ. Saf. 98, 162–170. doi: 10.1016/j.ecoenv.2013.09.007
Larson, C. A., and Passy, S. I. (2012). Taxonomic and functional composition of the algal benthos exhibits similar successional trends in response to nutrient supply and current velocity. Microbiol. Ecol. 80, 352–362. doi: 10.1111/j.1574-6941.2012.01302.x
Lecointe, C., Coste, M., and Prygiel, J. (1993). “OMNIDIA” software for taxonomy, calculation of diatom indices and inventories management. Hydrobiologia 269/270, 509–513. doi: 10.1007/BF00028048
Legendre, P., and Gallagher, E. D. (2001). Ecologically meaningful transformations for ordination of species data. Oecologia 129, 271–280. doi: 10.1007/s004420100716
Lengyel, E., Padisák, J., and Stenger-Kovács, C. (2015). Establishment of equilibrium states and effect of disturbances on benthic diatom assemblages of the Torna-stream, Hungary. Hydrobiologia 750, 43–56. doi: 10.1007/s10750-014-2065-4
Liu, J., Soininen, J., Han, B.-P., and Declerck, S. A. J. (2013). Effects of connectivity, dispersal directionality and functional traits on the metacommunity structure of river benthic diatoms. J. Biogeogr. 40, 2238–2248. doi: 10.1111/jbi.12160
Mangot, J.-F., Domaizon, I., Taib, N., Marouni, N., Duffaud, E., Bronner, G., et al. (2013). Short-term dynamics of diversity patterns: evidence of continual reassembly within lacustrine small eukaryotes. Environ. Microbiol. 15,1745–1758. doi: 10.1111/1462-2920.12065
Marcel, R., Bouchez, A., and Rimet, F. (2013). Influence of herbicide contamination on diversity and ecological guilds of river diatoms. Cryptogam. Algol. 34, 169–183. doi: 10.7872/crya.v34.iss2.2013.169
McArdle, B. H., and Anderson, M. J. (2001). Fitting multivariate models to community data: a comment on distance-based redundancy analysis. Ecology 82, 290–297. doi: 10.1890/0012-9658(2001)082[0290:FMMTCD]2.0.CO;2
Morin, S., Bottin, M., Mazella, N., Macary, F., Delmas, F., Winterton, P., et al. (2009). Linking diatom community structure to pesticide input as evaluated through a spatial contamination potential (Phytopixal): a case study in the Neste river system (South-West France). Aquat. Toxicol. 94, 28–39. doi: 10.1016/j.aquatox.2009.05.012
Needleman, S. B., and Wunsch, C. D. (1970). A general method applicable to the search for similarities in the amino acid sequence of two proteins. J. Mol. Biol. 48, 443–453. doi: 10.1127/1863-9135/2007/0168-0179
Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., et al. (2017). Vegan: Community Ecology package. R package version 2.4–3. doi: 10.1016/0022-2836(70)90057-4
Pandey, L. K., Bergey, E. A., Lyu, J., Park, J., Choi, S., Lee, H., et al. (2017). The use of diatoms in ecotoxicology and bioassessment: insights, advances and challenges. Water Res. 118, 39–58. doi: 10.1016/j.watres.2017.01.062
Passy, S. I. (2007). Diatom ecological guilds display distinct and predictable behavior along nutrient and disturbance gradients in running waters. Aquat. Bot. 86, 171–178. doi: 10.1016/j.aquabot.2006.09.018
Patrick, R. (1961). A study of the number and kinds of species found in rivers of the Eastern United States. Proc. Acad. Natl. Sci. Phila. 113, 215–258.
Peres-Neto, P. R., Legendre, P., Dray, S., and Borcard, D. (2006). Variation partitioning of species data matrices: estimation and comparison of fractions. Ecology 87, 2614–2625. doi: 10.1890/0012-9658(2006)87[2614:VPOSDM]2.0.CO;2
Pinckney, J. L., Hagenbuch, I. M., Long, R. A., and Lovell, C. R. (2013). Sublethal effects of the antibiotic tylosin on estuarine benthic microalgal communities. Mar. Pollut. Bull. 68, 8–12. doi: 10.1016/j.marpolbul.2013.01.006
Pont, D., Rocle, M., Valentini, A., Civade, R., Jean, P., Maire, A., et al. (2018). Environmental DNA reveals quantitative patterns of fish biodiversity in large rivers despite its downstream transportation. Sci. Rep. 8:10361. doi: 10.1038/s41598-018-28424-8
Proia, L., Lupini, G., Osorio, V., Pérez, S., Barceló, D., Schwartz, T., et al. (2013a). Response of biofilm bacterial communities to antibiotic pollutants in a Mediterranean river. Chemosphere 92, 1126–1135. doi: 10.1016/j.chemosphere.2013.01.063
Proia, L., Osorio, V., Soley, S., Köck-Schulmeyer, M., Pérez, S., Barceló, D., et al. (2013b). Effects of pesticides and pharmaceuticals on biofilms in a highly impacted river. Environ. Pollut. 178, 220–228. doi: 10.1016/j.envpol.2013.02.022
Ricciardi, F., Bonnineau, C., Faggiano, L., Geiszinger, A., Guasch, H., Lopez-Doval, J., et al. (2009). Is chemical contamination linked to the diversity of biological communities in rivers? Trends Anal. Chem. 28, 592–602. doi: 10.1016/j.trac.2009.02.007
Rimet, F., and Bouchez, A. (2011). Use of diatom life-forms and ecological guilds to assess pesticide contamination in rivers: lotic mesocosm approaches. Ecol. Indic. 11, 489–499. doi: 10.1016/j.ecolind.2010.07.004
Rimet, F., and Bouchez, A. (2012). Life-forms, cell-sizes and ecological guilds of diatoms in European rivers. Knowl. Manag. Aquat. Ecosyst. 406, 1–14. doi: 10.1051/kmae/2012018
Rimet, F., Chaumeil, P., Keck, F., Kermarrec, L., Vasselon, V., Kahlert, M., et al. (2016). R-Syst::diatom: an open-access and curated barcode database for diatoms and freshwater monitoring. Database 2016:baw016. doi: 10.1093/database/baw016
Roberts, D. W. (2007). labdsv: Ordination and Multivariate Analysis for Ecology. R package version 1.3–1.
Sabater, S. (2000). Diatom communities as indicators of environmental stress in the Guadiamar River, S-W. Spain, following a major mine tailings spill. J. Appl. Phycol. 12, 113–124. doi: 10.1023/A:1008197411815
Schloss, P. D., Westcott, S. L., Ryabin, T., Hall, J. R., Hartmann, M., Hollister, E. B., et al. (2009). Introducing mothur: open- source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 75, 7537–7541. doi: 10.1128/AEM.01541-09
Schneider, S., and Lindstrøm, E. A. (2009). Bioindication in Norwegian rivers using non-diatomaceous benthic algae: the acidification index periphyton (AIP). Ecol. Indic. 9, 1206–1211. doi: 10.1016/j.ecolind.2009.02.008
Schneider, S. C., and Lindstrøm, E.-A. (2011). The periphyton index of trophic status PIT: a new eutrophication metric based on non-diatomaceous benthic algae in Nordic rivers. Hydrobiologia 665, 143–155. doi: 10.1007/s10750-011-0614-7
Schwarzenbach, R. P., Escher, B. I., Fenner, K., Hofstetter, T. B., Johnson, C. A., von Gunten, U., et al. (2006). The challenge of micropollutants in aquatic systems. Science 313, 1072–1077. doi: 10.1126/science.1127291
Sipibel Report (2016). 2011–2015 Effluents Hospitaliers et Stations d’Épuration Urbaines: Caractérisation, Risques et Traitabilité – Synthèse des Résultats de Quatre Années de Suivi, d’Études et de Recherche sur le Site Pilote de Bellecombe (in French), 174.
Sládečková, A., Marvan, P., and Vymazal, J. (1983). “The utilization of periphyton in waterworks pre-treatment for nutrient removal from enriched influents,” in Periphyton of Freshwater Ecosystems. Developments in Hydrobiology, Vol. 17, ed. R. G. Wetzel (Dordrecht: Springer).
Stamatakis, A. (2014). RAxML Version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313. doi: 10.1093/bioinformatics/btu033
Stenger-Kovács, C., Lengyel, E., Crossetti, L. O., Üveges, V., and Padisák, J. (2013). Diatom ecological guilds as indicators of temporally changing stressors and disturbances in the small Torna-stream, Hungary. Ecol. Indic. 24, 138–147. doi: 10.1016/j.ecolind.2012.06.003
Stevenson, R., and Smol, J. (2003). “Use of algae in environmental assessments,” in Freshwater Algae in North America: Ecology and Classification, eds J. D. Wehr and R. G. Sheath (San Diego, CA: Academic Press), 775–804. doi: 10.1016/B978-012741550-5/50024-6
Stoof-Leichsenring, K. R., Epp, L. S., Trauth, M. H., and Tiedemann, R. (2012). Hidden diversity in diatoms of Kenyan Lake Naivasha: a genetic approach detects temporal variation. Mol. Ecol. 21, 1918–1930. doi: 10.1111/j.1365-294X.2011.05412.x
Tapolczai, K., Bouchez, A., Stenger-Kovács, C., Padisák, J., and Rimet, F. (2016). Trait-based ecological classify cations for benthic algae: review and perspectives. Hydrobiologia 776, 1–17. doi: 10.1007/s10750-016-2736-4
Tapolczai, K., Bouchez, A., Stenger-Kovács, C., Padisák, J., and Rimet, F. (2017). Taxonomy-or trait-based ecological assessment for tropical rivers? Case study on benthic diatoms in Mayotte island (France, Indian Ocean). Sci. Total Environ. 607, 1293–1303. doi: 10.1016/j.scitotenv.2017.07.093
Tornés, E., Mor, J. R., Mandaric, L., and Sabater, S. (2018). Diatom responses to sewage inputs and hydrological alteration in Mediterranean streams. Environ. Pollut. 238, 369–378. doi: 10.1016/j.envpol.2018.03.037
Vasselon, V., Bouchez, A., Rimet, F., Jacquet, S., Trobajo, R., Corniquel, M., et al. (2018). Avoiding quantification bias in metabarcoding: application of a cell biovolume correction factor in diatom molecular biomonitoring. Methods Ecol. Evol. 9, 1060–1069. doi: 10.1111/2041-210X.12960
Vasselon, V., Domaizon, I., Rimet, F., Kahlert, M., and Bouchez, A. (2017a). Application of high-throughput sequencing (HTS) metabarcoding to diatom biomonitoring: do DNA extraction methods matter? Freshw. Sci. 36, 162–177. doi: 10.1086/690649
Vasselon, V., Rimet, F., Tapolczai, K., and Bouchez, A. (2017b). Assessing ecological status with diatoms DNA metabarcoding: scaling-up on a WFD monitoring network (Mayotte Island, France). Ecol. Indic. 82, 1–12. doi: 10.1016/j.ecolind.2017.06.024
Verlicchi, P. (2018). “Final remarks and perspectives on the management and treatment of hospital effluents,” in Hospital Wastewaters – Characteristics, Management, Treatment and Environmental Risks. The Handbook of Environmental Chemistry, ed. P. Verlicchi (Heidelberg: Springer), 231–238.
Verlicchi, P., Al Aukidy, M., and Zambello, E. (2012). Occurrence of pharmaceutical compounds in urban wastewater: removal, mass load and environmental risk after a secondary treatment—a review. Sci. Total Environ. 429, 123–155. doi: 10.1016/j.scitotenv.2012.04.028
Verlicchi, P., Al Aukidy, M., and Zambello, E. (2015). What have we learned from worldwide experiences on the management and treatment of hospital effluent? — An overview and a discussion on perspectives. Sci. Total Environ. 514, 467–491. doi: 10.1016/j.scitotenv.2015.02.020
Villeneuve, A., Bouchez, A., and Montuelle, B. (2011). In situ interactions between the effects of season, current velocity and pollution on a river biofilm. Freshw. Biol. 56, 2245–2259. doi: 10.1111/j.1365-2427.2011.02649.x
Visco, J. A., Gentil, L. A. P., Cordonier, A., Esling, P., Pillet, L., and Pawlowski, J. (2015). Environmental monitoring: inferring the diatom index from next-generation sequencing data. Environ. Sci. Technol. 49, 7597–7605. doi: 10.1021/es506158m
Wang, Q., Garrity, G. M., Tiedje, J. M., and Cole, J. R. (2007). Naive Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl. Environ. Microbiol. 73, 5261–5267. doi: 10.1128/AEM.00062-07
Wiest, L., Chonova, T., Bergé, A., Baudot, R., Bessueille-Barbier, F., Ayouni-Derouiche, L., et al. (2018). Two-year survey of specific hospital wastewater treatment and its impact on pharmaceutical discharges. Environ. Sci. Pollut. Res. 25, 9207–9218. doi: 10.1007/s11356-017-9662-5
Windler, M., Leinweber, K., Bartulos, C. R., Philipp, B., and Kroth, P. G. (2015). Biofilm and capsule formation of the diatom Achnanthidium minutissimum are affected by a bacterium. J. Phycol. 51, 343–355. doi: 10.1111/jpy.12280
Keywords: pharmaceuticals, diatom communities, WWTP effluents, DNA metabarcoding, functional traits, water quality index, indicator species analysis
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
Received: 15 December 2018; Accepted: 15 March 2019;
Published: 09 April 2019.
Edited by:
Ghiglione Jean-Francois, Centre National de la Recherche Scientifique (CNRS), FranceReviewed by:
Gabor Varbiro, Centre for Ecological Research (MTA), HungaryYonghong Bi, Institute of Hydrobiology (CAS), China
Copyright © 2019 Chonova, Kurmayer, Rimet, Labanowski, Vasselon, Keck, Illmer and Bouchez. 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.
*Correspondence: Teofana Chonova, Teofana.chonova@gmail.com