Seasonal Variations and Resilience of Bacterial Communities in a Sewage Polluted Urban River

The Zenne River in Brussels (Belgium) and effluents of the two wastewater treatment plants (WWTPs) of Brussels were chosen to assess the impact of disturbance on bacterial community composition (BCC) of an urban river. Organic matters, nutrients load and oxygen concentration fluctuated highly along the river and over time because of WWTPs discharge. Tag pyrosequencing of bacterial 16S rRNA genes revealed the significant effect of seasonality on the richness, the bacterial diversity (Shannon index) and BCC. The major grouping: -winter/fall samples versus spring/summer samples- could be associated with fluctuations of in situ bacterial activities (dissolved and particulate organic carbon biodegradation associated with oxygen consumption and N transformation). BCC of the samples collected upstream from the WWTPs discharge were significantly different from BCC of downstream samples and WWTPs effluents, while no significant difference was found between BCC of WWTPs effluents and the downstream samples as revealed by ANOSIM. Analysis per season showed that allochthonous bacteria brought by WWTPs effluents triggered the changes in community composition, eventually followed by rapid post-disturbance return to the original composition as observed in April (resilience), whereas community composition remained altered after the perturbation by WWTPs effluents in the other seasons.


Introduction
Rivers flowing through cities are often used as receiving body for treated and untreated urban wastewaters all over the world [1,2]. Some of these sewage-contaminated rivers are amongst the most extreme examples of ecosystems disturbed by human activities. The impact of wastewater release on the receiving river depends primarily on the size of the city, the type of treatment applied to wastewaters and the flow of the river; the higher the river flow is, the higher its dilution capacity is. The UN predicted that, for 2050, 70% of the world population will be living in cities and this urban growth will mainly occur in less developed countries, where wastewater treatment facilities are scarce [2]. By 2025, it is expected that 27 megacities of more than 10 million inhabitants will exist, 21 of which in less developed countries [2]. Considering the importance, scarcity and fragility of freshwater ecosystems, research about the impact of the release of high amount of sewage on rivers functioning and ecological health is essential to preserve and have a rational management of this resource over a long term period.
It is now well known that microbial community inhabiting aquatic systems are one of the key players in the biogeochemical cycling of organic matter and nutrients and thus in the recovery and maintenance of ecosystem health and balance [3,4]. Sewage brings high loads of organic and inorganic pollutants that strongly modify the habitat of native freshwater bacteria. For instance, basic parameters such as temperature, pH and dissolved oxygen which are known as drivers of the bacterial community composition (BCC) [5,6,7] in aquatic systems can be influenced by urban wastewater discharges. Additionally, the loads of toxic trace metals from industries, of Persistent Organic Pollutants (POPs) and pharmaceuticals resulting from human activities could also affect the BCC [8,9,10,11]. Additionally, high concentrations of particulate matter within wastewater can modify light penetration in the water column and prevent phototrophic organisms' growth.
Sewage also brings high loads of allochthonous microorganisms [12]. If the allochthonous bacteria remain active after their release in the environments, as it has been previously suggested [13,14,15], they can have a major impact on the BCC of the river and could influence biogeochemical cycles. Besides, utilization of waters contaminated with allochthonous pathogenic microorganisms causes health risk.
Investigating the microbial population dynamic and composition of sewage polluted rivers and evaluating the resistance (insensitivity to disturbance) and resilience (the rate of recovery after disturbance (of microbial communities are essential for understanding how ecosystem health and functioning could be affected by the disturbances related to major sewage discharges. Despite the importance of these ecosystems and the increasing interest and advances in freshwater microbial ecology [16,17,18,19], papers devoted to BCC of urban rivers are still underrepresented in the scientific literature. Until now, the majority of the studies dealing with microbial community in urban aquatic settings have focused on indicator bacteria related to fecal and organic pollution [20,21,22,23]. Although some studies were done on polluted river sediments [24,25,26,27], to the best of our knowledge only few papers have examined the BCC of urban rivers waters [7,28,29,30]. The objectives of the present study were thus: (1) to study the spatial and temporal variations of the BCC along a sewage polluted urban river; (2) to identify the main environmental factors driving the microbial population dynamics in the watercourse; (3) to compare the river BCC before and after the discharge of major wastewater treatment plants (WWTPs) effluents; and (4) to analyze the fate and dynamics of different functional groups of sewage related bacteria once discharged into the river.
For this purpose, the Zenne River in Brussels (Belgium), a river highly disturbed by urban wastewater, was chosen as a model. This small river, with an annual average discharge of 4 m 3 s 21 upstream from Brussels, crosses the city and receives the effluents released by the two WWTPs of Brussels (treating together more than 1.5 million inhabitant equivalents); the average flow of treated wastewaters is in the same order of magnitude than the Zenne River flow rate upstream from Brussels. The river waters downstream from Brussels are thus roughly half composed of treated wastewaters, this proportion being even higher during the low flow periods of the river [31]. Seven stations along the river and the treated effluents of the two Brussels's WWTPs were sampled during four sampling campaigns at different seasons in 2010. The BCC was analyzed using 16S rRNA gene amplicon 454 pyrosequencing; in parallel, a large set of environmental parameters were measured.

Ethics Statement
No specific ethical or institutional permits were required to sample the Zenne River and the experimental studies did not involve endangered or protected species. We also obtained permission from Vivaqua (www.vivaqua.be) to sample the effluent of Brussels South WWTP and from Aquiris (www.aquiris.be) to sample the effluent of Brussels North WWTP.

Study Site and Sampling Strategy
The Zenne River is located in the Scheldt watershed and flows in the three Belgian regions (Wallonia, Flanders and Brussels-Capital); it is a tributary of the Dijle River (Fig. 1). The Zenne watershed (surface 991 km 2 ) is characterized by agricultural activities in its upstream part and an important urbanization in its downstream part. The annual average discharge of the Zenne River arriving in Brussels from the upstream part of the watershed is 4 m 3 s 21 (deduced from the last 20 years water flow data monitored by the Hydrologisch Informatie Centrum -HIC). Before the river reaches the Brussels area, it already receives several effluents of small scale WWTPs. The population density in the watershed is very high (on average 1260 inhabitants km 22 ) and mostly located in Brussels city and suburbs. The Zenne River has a total length of 103 km and crosses Brussels city from South to North over a distance of about 20 km. For sanitation purposes, the Zenne River was covered in the 19th century (1867-1871) in most of the stretch located in the Brussels region. In Brussels area, the Zenne River receives the sewage waters from two large WWTPs: the Brussels South WWTP (360,000 equivalent-inhabitants) in operation since the year 2000 and the Brussels North WWTP (1.2 million equivalent-inhabitants) in operation since 2007. The Zenne River also receives waters from two tributaries around the Brussels region the Zuunbeek, and the Woluwe (Fig. 1). Some small tributaries located in the Brussels area are diverted in the sewer collectors so that their waters reach Brussels WWTPs.
The two Brussels WWTPs function according to different technologies. The Brussels South WWTP treatment line includes a primary settling stage and a secondary biological treatment (activated sludge). At Brussels North WWTP, there are two treatment lines. The first one (biological line) includes a primary settling stage followed by a modern tertiary treatment technology (simultaneous removal of biodegradable organic carbon, nitrogen and phosphorus by an activated sludge process). The other treatment line (rain line) runs in parallel to the biological line when the discharge reaching the WWTP is too high in wet weather situations; this rain line uses only a coagulation-assisted high-rate lamellar settling process. On an annual basis, the volume treated in the biological line accounts for roughly 90% of the total volume reaching the WWTP.
Four sampling campaigns were performed in 2010 (January, April, July and October). Overall these sampling campaigns took place during marked dry weather periods; the discharge of the Zenne River just upstream from Brussels was always low during these campaigns: Comparatively, the April campaign was characterized by the driest conditions, both upstream and downstream Brussels city. Regarding October, there was actually a small rainfall occurring in Central Belgium during the sampling campaign but this did not even cause combined sewer overflow (CSO) discharge in the Zenne.
Seven stations were sampled along the Zenne River ( Fig. 1) in the stretch located downstream from the confluence with its major right-bank tributary the Sennette. Accordingly, a kilometric scale along the longitudinal course of the river was defined; it is arbitrarily set at zero at Lembeek (station Z1) and increases from upstream to downstream. Stations Z1 (0 km) and Z3 (13 km) are located upstream from Brussels. Station Z5 (20 km) is located just downstream from the Brussels South WWTP effluents release. Stations Z7 (30 km) and Z9 (34 km) are located upstream and downstream from the Brussels North WWTP, respectively. Stations Z10 (39 km) and Z12 (51 km) are significantly downstream from the Brussels conurbation area. Besides, in each of the sampling campaigns, the effluents of both Brussels WWTPs were sampled (average daily samples collected with refrigerated automatic samplers) and analyzed.

Physico-chemical Analysis and Bacterial Enumeration
Temperature, pH, and electrical conductivity (EC) were measured directly on site using a portable WTW 340 multiprobe. Dissolved oxygen was measured on the spot with a WTW oxi 323 field probe. Total alkalinity was determined by titration with HCl. Inorganic nitrogen (NH 4 , NO 2 and NO 3 ) and phosphorus (PO 4 ) concentrations were determined by automatic colorimetric methods using a QuAAtro (Seal, Analytical) segmented flow analyzer system. Suspended particulate matter (SPM) was estimated as the weight of material retained on Whatman GF/F glass fiber filter (diameter 4.7 cm, particle retention size 0.7 mm) per volume unit after drying the filter at 105uC. For dissolved (DOC) and particulate organic carbon (POC) analysis, samples were filtered on GF/F glass fiber filters previously combusted at 550uC. DOC was measured in the filtrate while the material retained on the filter was used to estimate POC. Glassware receiving the water samples for DOC analysis was muffled at 550uC for 4 h after cleaning; samples for DOC analysis were preserved with H 3 PO 4 (0.1% final concentration). DOC concentrations were measured with a total organic carbon analyzer in which inorganic carbon is eliminated by bubbling in the presence of phosphoric acid and organic carbon is oxidized by UV promoted persulfate oxidation and the produced CO 2 detected by infrared spectrometry (Hiper-TOC, Thermo). The GF/F filters with the particulate material were first acidified to remove inorganic carbon. POC was analyzed with an elemental analyzer (Flash EA 1112, Thermo) where organic material is oxidized by catalytic flash combustion at 1020uC and the produced CO 2 is measured by thermal conductivity. The DOC analyzer was calibrated using standard sucrose solutions, the POC analyzer was calibrated with acetanilide powder.
Bacterial abundance was determined by epifluorescence microscopy at 10006magnification, after DAPI staining according to the procedure proposed by Porter and Feig [32].

DNA Extraction
Microbial biomass was collected from water fraction and concentrate by filtration. An aliquot (from 0.8 L to 5.0 L) of each sample was filtered in triplicate directly through 0.2 mm pore-size, 47-mm-diameter polycarbonate filters (Millipore, MA). All filters were placed in 3 mL lysis buffer (50 mM Tris-HCl, pH 8.0, 40 mM EDTA, 750 mM sucrose) and stored at 220uC until use. DNA was extracted using a combination of enzymatic cell lysis and cetyltrimethyl ammonium bromide (CTAB) extraction protocol as previously described [33]. Dry DNA pellets were finally rehydrated in 100 mL of 10 mM Tris-HCl buffer (pH 7.4) and further purified using Min Elute reaction cleanup kit (QIAGEN). DNA concentration and purity were then determined using a Nanodrop ND-2000 UV-Vis spectrophotometer (Nanodrop, DE). Purified DNA extracts were stored at 220uC until use.

Bacterial 16S rRNA Gene Tag-encoded FLX-titanium Amplicon Pyrosequencing
Bacterial tag-encoded FLX gene amplicon pyrosequencing (bTEFAP) analysis was carried out by means of a Roche 454 FLX instrument with titanium reagents. Titanium procedures were performed at the Research and Testing Laboratory (Lubbock, TX, USA) based upon RTL protocols (www. researchandtesting.com) as previously described [34,35]. DNA samples were diluted to a final concentration of 20 ng/mL prior to bTEFAP. One-step PCR with a mixture of Hot Start and HotStar high-fidelity Taq polymerases was used. The PCR primers for FLX amplicon pyrosequencing were chosen to span the variable regions V1-V3 in the 16S rRNA gene: 27F (59-GAGTTT-GATCNTGGCTCAG-39) and 519R (59-GWNTTACNGC-GGCKGCTG-39) for Bacteria. These selected primers cover about 78% of publicly available 16S rDNA sequences of Bacteria (TestPrime tool at SILVA webpage http://www.arb-silva.de/ search/testprime/). All sequences generated in this study can be downloaded from NCBI Short Read Archive, accession numbers: (SRP021937).

16S rRNA Data Processing
The Research and Testing pipeline (Lubbock, TX, USA) performed the denoising were performed on each region. Briefly data were decompressed and sequencing errors were reduced by trimming flows (i.e., denoising) and sequences by applying the following criteria: amplicons shorter than 200 bp in length, -Afterwards, using USEARCH [36], reads were dereplicated meaning they were clustered together into groups such that each sequence was an exact match to portion of the seed sequence for the cluster. Then, the seed sequence from each cluster was sorted by abundance, largest cluster to smallest cluster, without any minimum size restrictions on the clusters. Afterwards, putative chimeras (checked by using USEARCH [36]) were removed from our data set. Then singletons were also removed. The denoised and chimera checked reads were condensed into a fasta file. Then failed sequences and sequences that have low quality tags, primers, or ends were removed.
The open-source, platform-independent, community-supported software program, Mothur (Mothur v.1.25.1; http://www. mothur.org) was used to analyze the clean sequence data. We trimmed the barcodes and primers from the resulting sequences.
Rarefaction curves based on identified OTU, Shannon diversity index and richness estimator Chao1 were generated using Mothur for each sample at 97% similarity. Sequences were classified using the greengenes database at 80% confidence threshold with Mothur. After taxonomic assignment of the sequences down to the phylum, class and genus level, relative abundance of a given taxonomic group was set as the number of sequences affiliated with that group divided by the total number of sequences per sample.

Statistical Analyses
Genera abundances were standardized, square root transformed, and assembled into a Bray Curtis similarity matrix to generate a multidimensional scaling (MDS) plot, using PRIMER 6 [37] and plotted to assess similarities among samples in twodimensional space. MDS plots represent relative distances among samples in relation to the rank order of their relative similarities. The Spearman correlation test was used to analyze the correlations between relative abundance of genera and phylum/ environmental parameters. The discrimination of bacterial assemblages based on water type was tested with one-way analysis of similarities (ANOSIM) using PRIMER 6 [37]. Similarity percentage (SIMPER) [38] was used to determine which genus contributed most to the dissimilarity between upstream stations (Z1, Z3), downstream stations (Z5, Z7, Z9, Z12 and WWTPs samples).
Rest of the analyses was conducted in R [39], using the vegan [40] and labdsv [41] packages. Procrustes analyses were done to determine the degree of concordance between MDS ordinations obtained at different taxonomic resolution levels (phylum, genera and OTU 0.03 ) on R [39]. Overall patterns in the microbiome were illustrated using heatmaps. When constructing the heatmaps, differences among sites were characterized using Bray-Curtis distances and differences among taxa were characterized using Euclidean distances. In addition, data within samples were standardized to facilitate comparisons. Species that occurred only once were removed for clustering analysis and only those having a relative abundance of more than 5% were shown in the heatmaps. We quantified the common and unique influences of environmental variables (T + pH + O2 + EC + SPM + NO2 + NO3 + NH4 + POC + DOC + PO4 + Alkalinity), distance and season on the bacterial community variation using variation partitioning analyses [42], using the varpart function in the R package 'vegan' package.

Physico-chemical and Biological Characterization of the Zenne River Samples
The four sampling campaigns were characterized by very contrasting temperatures ( Fig. 2A), covering a range from 4uC (January) to 22uC (July). An increase of the temperature is systematically observed after the discharge of WWTPs effluents regardless of the season. This effect is particularly marked after Brussels-North WWTP considering its high volumes of effluent discharge into the river.
In fall and winter, the amount of dissolved oxygen (DO, Fig. 2B) in water was quite high at the entrance of the study zone and decreased in the Brussels area. In winter a rapid decrease was observed after the discharge of the Brussels North WWTP, whereas in autumn the decrease was continuous all along the river. In contrast, in summer, DO was already very low at the entrance of the study zone (around 2 mg L 21 ) and dropped to 1 mg L 21 in the Brussels area values which are below the lethal concentration for fish [43,44]. The release of well aerated effluents from Brussels North WWTP obviously increased the oxygen level in the river that decreased rapidly further downstream. During our sampling campaigns, the particulate organic carbon (POC) content was detected in the range of 2-11 mg/L in all the stations, except three samples: Z5 in April, (51 mg/L, downstream from the Brussels South WWTP discharge) and Z10 and Z12 in January (23 mg/L, downstream from the Brussels North WWTP discharge). These abnormal increases in POC concentrations seem to be associated in April with the discharge of sludge from the activated sludge process at Brussels South WWTP. Such sludge release often occurs at Brussels South WWTP and was confirmed by microscopic observation for the studied April 2010 situation. In addition, in April, the driest studied period, the water mass at Z5 was very heterogeneous (due to the absence of a good mixing between WWTP South effluents and river water) so that particles were overrepresented in this sample. In January, the increase in POC downstream from Brussels area seems to be due to some untreated wastewater released into the river from station Z10. A significant increase in the ratio POC/SPM at Z10 and Z12 stations (25 to 30% C compared to 15% in Z9) was an indication of a high release of untreated sewage (with SPM rich in organic carbon) upstream from station Z10. On the other hand, DOC content was detected in the range of 5-13 mg/L in all the sampling stations. An increase in DOC concentration along the river was observed (Fig. 2D). Furthermore, in all the parameters measured, only NH 4 and PO 4 were positively correlated with each other, regardless of sampling season. These two parameters had a peak concentration downstream from the Brussels South WWTP (Z5) during each sampling campaign, indicating a major load of NH 4 and PO 4 from this WWTP effluent, due to the lack of tertiary treatment in Brussels South WWTP (Fig. 2E and 2F).
The total number of bacteria enumerated by epifluoresence microscopy ranged from 2.2610 9 to 2.3610 10 cells/L (data not shown). These bacterial abundances are in the upper range of the values reported in freshwaters but abundances of bacteria higher than 10 10 cells/L were already reported in rivers receiving high amounts of wastewaters [45]. Although ANOVA analysis revealed no significant difference between seasons, median values in January and April were higher than in July and October.  First of all, the data from river samples were compared to those from wastewater effluents, regardless of sampling time and station. Richness was estimated to be 1786 and 4407 OTUs 0.03 for WWTPs effluents and the river, respectively, while Chao indices were found to be more than 2000 for WWTPs samples, and 4500 for the river samples ( Table 1). The rarefaction curves at 0.03 (Fig.  S1A) showed that the global diversity of the two series of samples was well represented with the number of sequences analyzed. A recent study done on different treated wastewater effluents reported that species richness might vary between 1500 and 4000 based on the type of the treatment applied [46]. Detection of lower bounds for species richness of WWTPs effluents in comparison to river samples might be due the massive enrichment of some species which are very well adapted to the conditions encountered in some wastewater treatment processes such as activated sludge. When the samples were analyzed individually, richness was estimated between 573 and 1073 in January (Fig.  S1B). In April, samples reached to the saturation between 132 and 1271, increase in richness was observed after WWTPs discharges (Fig. S1C). In July, no sample except Z3 has reached to saturation (Fig. S1D). Richness was estimated between 349 and 774 in October (Fig. S1E). The highest richness in WWTPs effluents was detected in Brussels North in January (Fig. S1F). Figure 3 shows the taxonomic assignations at phylum/class level for all the 28 river samples together (A), and the 8 treated wastewaters samples (B). In the river (Fig. 3A), the dominant phylum was the Proteobacteria with up to 73% of the sequences, followed by Bacteroidetes up to 46%, Cyanobacteria up to 33%, Actinobacteria up to 21% and Firmicutes up to16%. Among the Proteobacteria, the most abundant class was the Betaproteobacteria (up to 63%), followed by the Epsilonproteobacteria (up to 26%) Gammaproteobacteria (up to 22%), and Alphaproteobacteria (up to 8%). Previous studies also reported the dominance of Proteobacteria in rivers [7,28] and the class Betaproteobacteria is known to be a typical freshwater clade [47]. At genus level (data not shown), 18% of the sequences remained unclassified. The most abundant genus was: Flavobacterium (Bacteriodetes) that represented 13% of the sequences, followed by Arcobacter and Malikia (both 7%, Epsilonand Betaproteobacteria, respectively) and Rhodoferax (4%, Betaproteobacteria). Sixteen percent of the total sequences corresponded to genera that represented less than 1% of the sequences from the river samples.
Concerning the treated wastewater effluents (Fig. 3B), at phylum level the global composition of the bacterial community was very similar to the one observed in the river samples. The dominant phylum was also the Proteobacteria with up to 74% of the sequences, followed by Bacteroidetes (up to 37%), Actinobacteria (up to 18%). Among the Proteobacteria, the most abundant class was the Betaproteobacteria (up to 43%), followed by the Alpha-and Epsilonproteobacteria (up to 16%), Gammaproteobacteria (up to 23%). Although community structures of river and wastewater effluents samples were pretty similar at phylum level, BCC of the treated effluents and the upstream samples were found significantly different at the genus level (R = 0.472, p,0.05). Therefore, after the release of the WWTPs effluent to the river, the BCC of the river might alter significantly. Recently, Drury et al. [48] studied BCC of rivers sediments and reported that WWTP effluents had impact on bacterial community structure and contributed to a biotic homogenization of rivers with different chemical and physical characteristics.

Multivariate Analysis of the Bacterial Community of All the River Samples
MDS analysis at OTU 0.03 and genus level was performed to see the BCC changes along the river at the different seasons (Fig. 4). Procrustes analyses showed that bacterial diversity patterns were strongly reproducible at both taxonomic levels (correlation coefficient = 0.75, p = 0.001) Despite the expected important influence of WWTPs discharge, the two major clusters based on the observed BCC corresponded to, on one hand, the spring and summer samples and, on the other hand, the autumn and winter samples. Previous studies on different river systems also showed large seasonal differences in BCC and bacterial production [49,50].
In addition, the samples from spring and summer are more scattered than those from autumn and winter suggesting more important variations along the course of the river in the warmer months.
The environmental parameters were shown to be significantly correlated with the samples distribution (Fig. 4A). Temperature was found to be the strongest parameter driving the BCC as already reported in previous studies [50]. Other important parameters were: oxygen concentration, DOC, the SPM content and inorganic nitrogen (NH 4 , NO 2 and NO 3 ) concentrations. These results suggest that the marked seasonality of bacterial community could be associated with in situ bacterial activities (DOC and POC biodegradation associated with oxygen consumption and N transformation) which are higher in spring and summer. Indeed, BCC in April and July were positively correlated with PO 4 , NH 4 , total alkalinity and negatively with oxygen , showing active associated respiration processes, whereas opposite interactions were demonstrated in January and October. In situ bacterial activities in the river are lower in the colder situations; therefore the input of allochthonous bacteria by WWTPs effluents had a stronger impact on the BCC of the river in fall and winter.
Furthermore, clear association between the phylotypes and the samples were shown (Fig. 4B). An important number of taxa were associated with autumn and winter samples while a lower number were associated with spring and summer samples. Four out of the five genera which were significantly correlated with the spring and summer samples: Planktoluna, Limnohabitans, Flavobacterium, Anabaena, have been clearly identified as planktonic bacteria autochthonous from freshwaters [51]. Among the numerous genera present in winter and autumn samples, several belong to groups associated with sludge and wastewaters (Faecalibacterium, Haliscomenobacter, Beijerinckia, Rhodocyclus, Prevotella, Roseburia, Nostocoida).
This supports the hypothesis of a higher importance of in situ microbial processes in spring and summer, carried by taxa well adapted to the environmental conditions in the river, while allochthonous bacteria associated with wastewaters were found dominant during winter and autumn.
To estimate the relative importance of environmental, seasonal and spatial components in shaping bacterial community structure, we used a variation partitioning analysis. The percentages of the total variation in the bacterial matrix that can be attributed to the different components of variation (environmental, spatial, seasonal and unexplained) were based on the adjusted (unbiased) fractions. When considering the total bacterial diversity, the highest part of the variation (65%) was remained unexplained. The highest part of variation was explained by environmental factors (29%). Despite the important discharges of wastewaters along the river and the significant modification of the river BCC due to the allochthonous bacteria inputs, followed by the strict influence of the sampling location and season (3% each); Besides, correlation between OTUs 0.03 and environmental parameters were calculated (rho.0.8, p,0.05), and POC, COD and SPM showed significant correlations with some species. For instance, Acinetobacter, Rhodococcus corprophilus, Thauera were found to have a linear relationship with POC (Table 2). Such correlations were not previously shown at a lower taxonomic level. However, a significant correlation between Bacteriodetes and POC was also previously shown in a study conducted in the English Channel [52].

Spatio-temporal Variations of Bacterial Community Structure
Seasonal fluctuations of bacterial community. In this section, all the samples from each sampling campaign have been grouped and the four groups have been compared in order to evaluate the seasonal impact. Figure 5 presents box plots of Chao and Shannon indices for the different seasons (all sampling stations of each campaign grouped). The median of both parameters was the lowest in April; and a gradually increasing trend was observed from April to January (Fig. 5A). The richest and most diverse samples were those from January. Besides, in January the Shannon index had very low standard deviation, while the highest deviation was observed in April (Fig. 5B).
The medians of Chao indices ranged from 600 (April) to 1300 (January), it means that more than twice different species are expected to be found in January with regards to April. The Shannon index also varied a lot ranging from 3.4 (April) to 5.2 (January). There are different seasonal studies done in the aquatic environments showing higher richness estimations in winter [53,54]. Due to the lower bacterial activity in colder situations, species are more evenly distributed in colder season than in the warmer seasons. However, temperature cannot be the only parameter explaining the difference in diversity indices. Combi-nation of different parameters such as oxygen, organic carbon and/or nitrogen components might have influenced the richness.
Accordingly, effect of seasonal variation was also observed on BCC. ANOSIM analysis (Table 3) revealed that the BCC of the samples from each season were significantly different from each other with the exception of BCCs in April and July which did not differ significantly. Besides, the SIMPER analysis was used to determine the percentage of similarity between the four sampling campaigns in the upstream stations (Z1 and Z3), in order to test if there was any seasonal impact on the BCC of upstream waters before receiving Brussels WWTPs effluents. SIMPER analysis also revealed which individuals contributed most to the dissimilarity between samples (Table S1). The data showed that BCC of the upstream samples in cold seasons (January and October) had the highest similarity with 62%, followed by similarity between April and July with 53%. On the contrary, the highest difference in the BCC was observed between January and April with 63%. Bloom of Malikia and Flavobacterium in April and relatively high abundant taxa such as Dechloromonas, Lepthothrix, Pseudomonas and Arcobacter in January explained the difference between the BBC of the two seasons.
Impact of WWTPs effluents on the river BCC. As already mentioned, WWTPs effluents have not the same characteristics in terms of BCC than the receiving water body; therefore, it is important to evaluate the impact of WWTPs effluents on the BCC of the river. First of all, we analyzed if the season has an impact on the BCC of the WWTPs' effluents. Based on the richness and diversity analysis, Brussels South WWTP had clearly the lowest richness and diversity index in July. Besides, BCC analysis revealed that Brussels South WWTP in July clustered apart from other WWTP effluents (Fig. S2). ANOSIM revealed no further clear separation between the WWTPs effluents based on the season or WWTP (p.0.05). Figure S3 presents the longitudinal profiles of the Chao and Shannon indexes, each point representing the median of the four sampling seasons. The richness increased after the release of Brussels North WWTP, the standard deviations were large for all the sampling stations, reflecting the important seasonal variations. The diversity (Shannon index) slightly increased after the release of both WWTPs and these stations located just downstream from the WWTPs presented lower standard deviations (indicating lower seasonal variations) than the other sampling stations. To evaluate the impact of WWTPs discharge, three artificial groups of samples were compared: (1) the 8 samples of treated wastewaters; (2) the 8 river samples collected upstream from Brussels South WWTP discharge (Z1 and Z3 stations from the four sampling campaigns); (3) the 20 river samples collected downstream from Brussels South WWTP discharge (5 samples per sampling campaign). First of all, the richness estimation did not change significantly after the discharge, although significant effect of WWTPs discharges on the BCC in the water column was observed. In contrast, a recent study  showed that WWTP effluent significantly reduced bacterial diversity and richness in rivers sediments [48], which might indicate a different response of BCC to the WWTP discharges in the water column and the sediments. It has also been previously reported that the community structure of free living bacteria and particle-attached bacteria were different [16]. Due to the settling of the bacteria attached to particulates, the impact of WWTPs effluents on the bacterial richness in the water column and the sediments might be different. Besides, ANOSIM analysis (Table 4) revealed that upstream samples were significantly different from downstream and WWTP effluents, while no significant differences were found between WWTPs effluents and the downstream samples. SIMPER analysis (Table S2) shows that no single genus can explain more than 5% of the differences. This means that the differences between the upstream and downstream samples are due to a shift after the discharge of very diverse bacterial community from the WWTPs. The impact of WWTPs on the BCC was studied further per season.
Longitudinal profiles of BCC in the river at different seasons. As January and April presented the most contrasting seasonal impact on the bacterial community, they are discussed in more detail below. Longitudinal profiles of BCC in July and October are further discussed in Results S1 and their dendrograms were illustrated in Figure S4.
In January, two major clusters were observed: river samples (Z1 to Z9), on one hand, and wastewater effluents and river sample Z12, in which a release of raw wastewater is suspected, on the other hand (Fig. 6A). River samples had three subclusters: samples collected upstream of the Brussels South WWTP (Z1, Z3); samples  collected downstream from the Brussels South WWTP (Z5, Z7) and samples collected downstream from the Brussels North WWTP (Z9, Z10). It is interesting to note that the BCC of the samples taken just downstream from the WWTPs effluents discharge were not dominated immediately by the taxa brought by the effluents. The higher abundance of Intrasporangiaceae (closest match to Tetrasphaera), Saprospiraceae (closest match to Haliscomenobacter), and Rhodoferax originating from Brussels North WWTPs effluents were only observed at Z12 station. Gram-positive Tetrasphaera-related organisms (Actinobacteria) are putative polyphosphate-accumulating organisms that are abundant in many fullscale enhanced biological phosphorus removal plants [55] and Haliscomenobacter are one of the filamentous bacteria detected in wastewater treatment plants which are responsible for biomass bulking [56,57]. Their survival in the cold temperature was previously reported [57,58].
In contrast to January, a strong resilience of bacterial communities in the river was observed in April. Bacterial community was strongly modified just downstream from the WWTPs effluents discharges (Fig. 7) and in few kilometers the initial bacterial community composition was recovered and this phenomenon was observed after both WWTPs discharges. This difference between BCC in April and January was also shown by clustering. Two major clusters were observed in April, on one hand, the WWTPs samples and the river samples collected just downstream from the effluents and, on the other hand, the remaining river samples (Fig. 6B). The system seems to be resilient, since the initial situation is recovered already a few kilometers downstream from the discharges of WWTPs (disturbance). SIMPER analysis revealed the dissimilarity between WWTPs and upstream samples 71%. This was due to the bloom of Malikia and Flavobacterium in the upstream stations and the dominance of Zooglea, Arcobacter, Acidovorax, and Acinetobacter in the WWTPs' effluents. Arcobacter and Zooglea were previously isolated from sewage [59,60]. Flavobacteria are found in different environment such as soil and freshwater, some of whose species are known to cause disease in freshwater fishes [61] whereas hydrogen-oxidizing and polyhydroxyalkanoate accumulating Malikia are restricted to freshwater habitats [62]. We could not find any correlation between Malikia bloom and any environmental parameters measured in April. Malikia bloom might be due to the low level of river flow at this driest sampling time. Besides, in contrast to January, the impact of WWTPs effluents release has only been observable just downstream from the WWTPs discharges (Z5 and Z9 stations), while the BCC of samples collected further downstream from the WWTPs discharges were similar to the BCC in the upstream part of the river (Fig. 6B). In this case, the composition and structure of bacterial community seems to be dominated by in situ processes carried by bacterial groups well adapted to the riverine conditions, with an exception of the sampling points just downstream from the WWTPs discharges.
As it is described above, the impact of WWTPs effluents was season dependent. A recent study showed that sediments of two different rivers with different characteristics in their upstream watershed had almost indistinguishable BCC downstream from the WWTPs discharges, suggesting the potential of WWTPs to reduce the natural variability [48]. Our data suggest that the resilience of the BCC in the water column is season dependent. Studies done by Yannarell et al. [63] and Crump et al. [50] had also contrasting results based on the seasonal impacts on the variability of bacterial communities, which indicates the influence of environmental factors other than temperature like in our case.

Conclusions
Our data demonstrated that BCC was directly or indirectly influenced by complex seasonal shifts in environmental parameters. Our study showed that disturbance strength has the potential to change the composition of bacterial communities and confirms previous findings that bacterial communities are generally not resistant to disturbances [48,64,65]. Indeed, community composition remained altered after the perturbation by WWTPs discharges in three out of the four seasons. However in April,the driest sampled season -, recovery to pre-disturbance community composition was observed just few kilometers downstream from the WWTPs (pulse disturbance). This is an evidence of the bacterial community resilience in a strongly polluted urban river. Studies dealing with in-situ evidences of freshwater bacterial community stability are scarce and essential for accurate modeling and management of water resources. Environmental and biological factors associated with resilience and resistance must be further studied to assess and prevent the impact of (treated or untreated) sewage pollution on the freshwater ecosystem functioning. When constructing the heatmaps, differences among sites where characterized using Bray-Curtis distances and differences among taxa were characterized using Euclidean distances. In addition, data within samples were standardized to facilitate comparisons. The relative values for bacterial genera are depicted by color intensity. (PDF)  Results S1 Longitudinal profile of BCC in the river in July and October.