Water Purifying Capacity of Natural Riverine Wetlands in Relation to Their Ecological Quality

Wetlands are among the crucial nature based solutions for river water purification. The ecological quality of these multi-purposed ecosystems could be disrupted due to overloading of nutrients and other pollutants. We investigated the water purifying efficiency and ecological quality of three natural riverine wetlands located in urban area in southern Ethiopia. Of the three wetlands, Boye wetland is located in a densely populated area with intensive human activity, and the other two, Fisho and Kitto located in a relatively less populated area of Jimma town. We sampled water, macroinvertebrates, and diatoms, to compare sites before joining the wetland, within the wetland and after passing through the wetland. Considering both seasons, up to 74% TP, 73% DIN, and 77% BOD reduction was recorded. The lower concentration of nutrients, and BOD in sites after joining the wetland showed the presence of pollution attenuation. Macroinvertebrate and diatom based bio-indices depicted higher biodiversity and lower relative abundance of tolerant taxa in sites after joining the wetland, which implied the potential of studied wetlands to reduce pollutants and sustain biodiversity. The incoming streams of Boye showed bad water quality and heavily degraded ecological status (Ethbios score 8–10). Most of the sites associated with Boye wetland depicted a major ecological disturbance (Ethbios score < 44). The incoming streams of Fisho revealed moderate (Ethbios score = 58) and poor (Ethbios score = 20) water quality. Most of Fisho sites had a moderate water quality with significant ecological disturbance. Sites associated with Kitto wetland, had a slight ecological disturbance (Ethbios score of 93 and 72) and some with significant ecological disturbance (Ethbios index score between 45 and 59). The wetland fed by heavily polluted streams showed the most degraded ecological quality compared to the other two that are fed by less polluted streams. An RDA model visualized pollution gradient among sites. Generally, this study confirmed the potential of natural wetlands to mitigate nutrients and organic pollutants and sustain biodiversity. However, when the incoming water is heavily degraded, the retention of pollutants seriously affect the wetland’s ecological quality.


INTRODUCTION
Freshwater pollution is one of the worsening environmental issues in the developing countries (Hanasaki et al., 2013;Babayemi et al., 2016), commonly associated with nutrient, organic and inorganic loadings, alteration of hydro-morphology, and habitat degradation (Chen et al., 2016;Sharifinia et al., 2016). Intensive agricultural practices with the application of uncontrolled fertilizers and pesticides and human wastes are the most important sources of freshwater nutrient loading, which mostly occurs in a diffused form and threaten the aquatic ecosystem (Strokal et al., 2016(Strokal et al., , 2017. Such non-point source pollutants are more challenging for both containment and treatment (Verhoeven et al., 2006). Hence, a sustainable water resource management which can uphold the societal needs while keeping ecological, environmental and hydrological integrity is crucial (Berhanu and Poulton, 2014;Serpa et al., 2014).
However, wetlands could be overloaded due to excessive nutrients and other contaminant inputs, thus the ecosystem may shift into another stable state; from being a sink to a source of pollutants (Verhoeven et al., 2006;Resende et al., 2010). This might result in loss of ecological quality and a shift in diversity, and composition of taxa (Hansson et al., 2005;Davidson et al., 2012).
In Ethiopia, wetlands are estimated to cover about 1.4% of the country's landmass (Gebresellassie et al., 2014). Like in many other African countries (Beuel et al., 2016), wetlands in Ethiopia are under continuous pressure of degradation (Mereta et al., 2012;Moges et al., 2016). Such problems are aggravated by rapid population growth in the absence of adequate sanitation facility and improper wastewater management (Bateganya et al., 2015). For instance, in eastern Ethiopia, people kept putting pressure on lake Haramaya until it vanishes (Mengisteab, 2012;Abebe et al., 2014). Other research reports indicated the degradation of, for instance "Boye" and Gilgel-Gibe wetlands (Van der Bruggen et al., 2009) and the complete dry out of "Aba Samueal" marsh (Gizaw et al., 2004). Other researchers also warn similar incidents could happen to other wetlands in the country if proper measures are not in place (Bekele et al., 2001;Beyene et al., 2012). An integrated approach should be used to assess ecological status and ensure the functioning of wetlands and, therefore to monitor and protect these multi-purpose ecosystems (Justus et al., 2016). In order to monitor the water quality, physicochemical parameters are important variables. They provide information about the status of the water during the sampling event (Resende et al., 2010). Additionally, biological indicators provide information on the general ecosystem and they can trace back past environmental variability. The biological integrity of an ecosystem is informative about the system's ability to support and sustain biodiversity (Lakew and Moog, 2015). Biological communities, particularly macroinvertebrates and diatoms, are effective to assess the level of disturbance and the overall ecological status of an aquatic ecosystem.
Wherever possible in water quality assessment, simultaneous use of two or more biological quality elements including macroinvertebrates, fishes, macrophytes, planktons, and diatoms is recommended to provide an integrated measurement of the overall ecological status. Researchers recommended the use of two or more suitable groups contemporaneously so that one can compensate the limitations of one group with the strength of the other (Carlisle et al., 2008;Torrisi et al., 2010;Lainé et al., 2014). However, for extremely high turbidity in the water system which is the case in our study area, fish and macrophytes are either absent or contain a very limited number of species to be used as bioindicators for ecological quality assessment (Krstic et al., 2007;Beyene, 2010). The simultaneous use of diatoms and macroinvertebrates reported to give good indication of ecological quality in freshwater (Resende et al., 2010;Justus et al., 2016).
Macroinvertebrates are indicators of choice to assess the presence of catchment disturbance and hydro-morphological alterations (Birk et al., 2012;De Troyer et al., 2016). Their ecology is well determined and responds sensitively to environmental gradients (Lakew and Moog, 2015). Diatoms community composition is more related to the water quality (Sonneman et al., 2001;Rimet et al., 2005). Different diatom species have their own distinct habitat preferences (Vilmi et al., 2015) which made them ideal to predict ecological quality. Beyene et al. (2009) suggested diatoms as powerful indicators of gradients of pollution in tropical freshwater ecosystems. These two bioindicators, macroinvertebrates and diatoms, are considered as complementary tools to assess water quality, and ecological status of the aquatic ecosystem (Feio et al., 2007;Resende et al., 2010).
This study used both bioindicators, diatoms, and macroinvertebrates, in relation with physicochemical variability to provide a complimentary insight on the ecological quality of the system under study. Since there is variation in hydrological variables during the wet and dry season, research in tropical water bodies should deal with seasonality (Beyene et al., 2014). In this study, we also consider the seasonal variation in water quality.
Although there are studies dealing with the water purification potential of wetlands (Vymazal, 2007;Fountoulakis et al., 2009;Mereta et al., 2012), both natural and constructed, or its ecological status, both of these aspects are seldom studied together. These two aspects should be studied together in order to design a wetland management tool which enhances the application of wetland ecosystems for in situ water purification without risking its ecological quality. Therefore, the main objective of this study is to determine the efficiency of natural riverine wetlands to reduce nutrient and organic pollutant concentrations, and to assess their ecological quality in relation to the incoming river water quality status. The specific objectives of the study are hence (a) to evaluate the efficiency of the three riverine wetlands in reducing pollution, (b) to assess the ecological status of incoming and outgoing streams to and from three natural riverine wetlands, and (c) to investigate the impact of the different incoming streams on the ecological quality of the wetlands themselves using physicochemical parameters and two (diatom and macroinvertebrates) biological quality elements.

Study Sites and Samples
The study site is located in Jimma, south western part of Ethiopia, between the geographic coordinates of 7 • 41 N and 36 • 50 E. The study area has an elevation of 1,704 m above sea level (Duguma et al., 2012). For this study, three riverine wetlands were selected namely; Boye (B), Kitto (K), and Fisho (F) wetlands (Figure 1). These wetlands are considered as urban influenced (Moges et al., 2016). However, the level of disturbance differs among the wetlands. Boye wetland is located in the most populated area of Jimma town with heavy anthropogenic activities. It is estimated to cover 60 ha area (Moges et al., 2016) on average. The major anthropogenic disturbances in this wetland, as observed during the field survey, include solid and liquid waste disposal, car washing, livestock grazing, and papyrus harvesting. There are four main streams flowing into Boye wetland, namely Kochie (B1), Dololo (B2), University (B3), and Awaitu streams (B4).
All of these streams serve as a waste dumping sites for the respective residents. They discharge all kinds of wastes without any form of prior treatment. These streams confluence within the wetland and flow to Gilgel Gibe River. Fisho and Kitto wetlands are located in a relatively less populated area of the town with different types of disturbances. The streams feeding them are also relatively less stressed compared to those that feed Boye wetland. Fisho wetland is estimated to cover 40 ha area. Two streams join Fisho wetland namely, Degoye 1 (F1) and Degoye 2 (F2). Brick mining, farming, clothes washing, and livestock grazing are the major disturbances in Fisho wetland. Kitto wetland with 50 ha estimated area, has two streams flowing into it. One is Kitto River (K1), located in the upstream of the wetland with most of its catchment upstream of the town, which is the least disturbed site, and the second stream (K2) is a merge of small streams passing through residential areas. The major anthropogenic activities in Kitto wetland site are cloth washing, swimming, and plant harvesting (especially papyrus). All the feeding streams were permanently flowing. However, their flow rate extremely decreased during the dry season.
A total of 20 sampling sites, including the streams associated with the wetlands, were selected to analyze the water purifying effect of the wetlands and their ecological status. Sampling locations before joining the wetlands were; B1, B2, B3, B4 (all four streams that join Boye wetland), F1, F2 (the two streams that feed Fisho wetland), K1 and K2 (the two streams feeding Kito wetland). Those with in the wetland were B5, B6, B7 (from Boye wetland), F3, F4 (from Fisho wetland), K3 and K4 (from Kito wetland). The sites immediately after passing through each FIGURE 1 | Map of the study area including sampling locations in three natural riverine wetlands of Jimma, southwest Ethiopia (Namely; Boye (Sampling points starting with "B"), kito (that start with "K"), and fisho (points starting with "F"); and the respective incoming streams that fed those wetlands together with the outflows of the three wetlands (see Supplementary Table S1 for the site code descriptions). The indicated flow direction represents all the three wetlands in the study map.
Frontiers in Environmental Science | www.frontiersin.org wetland were B8, F5, F6, K5, and K6. To deal with the seasonal variations in hydrological processes of tropical water bodies (Beyene et al., 2009), water samples, macroinvertebrates, and diatoms were sampled in replicates and duplicates, when access is difficult, in August (2016) and February (2017) representing the wet and dry seasons, respectively.

Measurement of Environmental Parameters
Dissolved oxygen (DO), water temperature, pH, and electrical conductivity were measured using a multi-probe meter (Hach-Model-HQ30d Single-Input, Multi-Parameter Digital Meter) and turbidity was measured using a Wag-WT3020 turbidity meter on-site. We filtered the water samples onsite using Whatman glass microfiber filters (25 mm Ø) and collected in polyethylene bottles of 100 ml for nutrient analysis (nitrite, nitrate, ammonium, and soluble reactive phosphorus). Wellmixed unfiltered water samples were collected using polyethylene bottles for total nitrogen (TN), total phosphorus (TP), total suspended and dissolved solids (TSS and TDS), chloride and alkalinity. The water samples for BOD analysis were collected by inserting clean polyethylene bottles of 1 L to a 30 cm depth in the opposite direction of the current flow of the rivers and immediately sealed and closed. The water samples were transported in a cold box to the laboratory of Environmental Health Sciences and Technology, Jimma University, Ethiopia. The samples were immediately put in to cold box and kept in freezer until the analysis. Generally, the samples were analyzed following standard methods described by APHA (2005). Environmental variables, total nitrogen (TN) and phosphorous [both soluble reactive (SRP) and total (TP)] were measured by Kjeldahl Nessler and ascorbic acid methods, respectively following APHA (2005). Biochemical oxygen demand (BOD) was measured using the azide modification of the Winkler's titrimetric method by determining DO contents of the samples before and after 5 days' incubation at 20 • C. High (0.3-30 mg/l) and low (0.01-0.5 mg/l) ranges of nitrate were measured using the cadmium reduction method of 8039 and 8192, respectively. We also used Nessler method 8038 for the determination of Ammonium-Nitrogen following the HATCH water analysis hand book of the USEPA available online for free 1 .

Macroinvertebrate Sampling and Identification
Macroinvertebrate habitats were detached by kicking the bottom sediment of the wetlands and in between vegetation and the river bed continuously for about 5 min within 10 m stretch. The dislodged macroinvertebrates were collected using a rectangular frame, 50 × 33 cm, kick net with a 250 µm mesh size. Then the content of the net was transferred into a white tray, sorted on-site using forceps and kept in vials containing 80% ethanol, and labeled with the site code. This procedure was performed twice for each of the 20 sampling locations, and samples were pooled together in a vial for each site. All the Diatom Sampling, Processing, and Identification Diatom assemblages were sampled simultaneously with macroinvertebrates and water samples. From each of the sampling locations two types of substrates, stone, and hardwood, were scraped using a toothbrush and pooled to a 50 ml polyethylene vial. In the absence of stones, we collected submerged vegetation to clean off the diatoms attached to it. Diatom samples were fixed with 4% formalin and transported to the biological laboratory of Jimma University, Department of Environmental Health Science and Technology.
Diatom samples for species level identification are mostly digested using concentrated H 2 SO 4 and HNO 3 according to Lange-Bertalot (2000) and Round et al. (2007). A concentrated acid oxidation method modified from Stevenson and Bahls (1999) was used for the cleaning of diatom frustules. The major modification in the diatom sample processing is the use of HCl as a pre-cleaning in order to remove calcium from the sample.
Slides were prepared by mounting 30 µl from a known volume of digested sample using Zyrax mounting media. Diatom identification and counting was done using a Carl Zeiss light microscope with 1,000× magnification. A minimum of 300 diatom valves were counted in each slide except for some samples from the wet season, where less than 300 diatom assemblages were found. The diatom valves were identified to the smallest possible taxonomic level using different standard identification keys (Lange-Bertalot, 1980, 2001, 2011Lange-Bertalot, 1986, 2000;Taylor et al., 2007).

Statistical Analysis
All the analyses were done using R statistical software (R Core Team, 2015, version 3.2.2) except for diatom indices, which was done using OMNIDIA version 5.2. Packages used in the analysis include Vegan (Oksanen et al., 2015), Packfor (Dray et al., 2013), Hmisc (Harrell et al., 2015).
Descriptive statistics were used to determine the overall status of the water chemistry along the sampling locations. All environmental parameters, except pH, were log-transformed [log (x + 1)] prior to analysis. Distribution of environmental parameters in each group was tested using Shapiro-Wilk test, which is the most powerful normality test (Razali and Wah, 2011). The normality was not attained in all groups, for the parameters under study, therefore non-parametric tests were used for the proceeding analysis. A Kruskal-Wallis rank sum test and Wilcoxon rank-sum test were used to determine the presence of a significant difference in environmental parameters among the three wetland sites and between the two seasons, respectively. The Kruskal-Wallis test was followed by Mann-Whitney test to make a pairwise comparison. Bonferroni correction was applied to correct for multiple testing in order to reduce the chance of type I error. Spearman correlation analysis was done to determine the association between environmental variables and biological indices using the Hmisc package.

Percent Change Calculation
The potential of wetlands to reduce or remove nutrients and organic pollutants was calculated, for TP, PO 4 -P, NO 3 -N, NH 4 -N, DIN, TN, EC, and BOD, using the following equation: where C in represent the incoming concentration and C out for the outgoing concentration.
For the wetlands having multiple inflows, the mean concentration of the inflows was taken as the incoming concentration (C in ).

Ordination Analysis
Macroinvertebrate taxa and diatom taxa abundance data were Hellinger transformed prior to proceeding analysis. In order to choose between linear, Principal Components Analysis (PCA), Redundancy Analysis (RDA) or unimodal; Correspondence Analysis (CA) models, Detrended Correspondence Analysis (DCA) was applied on the macroinvertebrate and diatom taxa count data. The maximum gradient length of the four ordination axis was between three and four in both cases (macroinvertebrates and diatoms), during both seasons (wet and dry season). The gradient didn't show the dominance of either linear or unimodal gradients. Hence, we applied both RDA and CCA to select the one which explains better, and RDA was chosen as it was the best from the two. RDA model was built by a forward selection procedure (permutation, n = 999) using double stopping criterion to investigate the relation between macroinvertebrate and diatom taxa with environmental variables. For each season, dry and wet, the RDA model was built for the three wetland sites together.

Compositional and Diversity Metrics
Macroinvertebrate and diatom-based compositional and diversity metrics were calculated to describe the community structure of diatoms and macroinvertebrates in each site, for both seasons.
Total taxa richness, Shannon's diversity, % filter collectors (FC), and Ethbios index were calculated based on macroinvertebrate community assemblage. Total taxa richness and %FC are among the recommended indices to assess the ecological quality of wetlands in Ethiopia (Mereta et al., 2013). Shannon's diversity is among the commonly used indices in ecological studies and it accounts for abundance and evenness (Magurran, 2013). Ethbios index is a recently developed index similar to biological monitoring working party (BMWP) but includes those macroinvertebrate taxa which are endemic to Ethiopia. The Ethbios is designed mainly for river systems of Ethiopia. To use it for the sampling sites from the wetlands we tried to verify it with the MMI developed for Ethiopian wetlands by Mereta et al. (2013). In determining the Ethbios index, each taxon is assigned a score based on its sensitivity to pollution; highest value for highly sensitive taxa and low values for pollution tolerant ones. The overall Ethbios index value of the site was obtained by summing up the scores for each taxon (Lakew and Moog, 2015).
Diatom based metrics were also used to have a comprehensive insight about the ecological quality of the system under study. The calculated indices include Shannon's diversity, total taxa richness, and the relative abundance of pollution tolerant taxa (% PT).

Purification Capacity of Natural Riverine Wetlands
The sites from the three wetlands were pooled together as those before joining the wetland, within the wetland, and after passing through the wetland. The pooled values of environmental parameters indicated the presence of pollution gradient. In the dry season, significantly lower concentration of TN (W = 0, p = 0.03), DIN (W = 1, p = 0.02), TP (W = 0, p = 0.03), and PO 4 -P (W = 0, p = 0.03), was observed in sites after draining through the wetland than sites before joining the wetland (Mann-Whitney Wilcoxon). Considerably higher concentrations of NO 3 -N and DO, and lower BOD was recorded within the wetland (see Supplementary Figure S1).
Percent reduction analysis exhibit variation among studied wetlands. During the dry season, Boye wetland showed the highest percent reduction of TP (74%), whereas Fisho and Kitto wetlands showed 23 and 26% reduction, respectively. The highest reduction of TN (43%) and BOD (18%) was recorded in Fisho wetland (Figure 2 dry season and Supplementary Table S3).
Concerning NH 4 -N and NO 3 -N, the highest reduction was 34 and 88% in Fisho and Kitto wetland respectively during the dry season (Supplementary Table S3). In the wet season, the highest percent reduction of BOD was recorded in Boye wetland; All results are in mean and SE (indicated in brackets) of all the sites in a wetland system per season. *Shows significant seasonal difference where, *P < 0.05, **P < 0.005, ***P < 0.001 (Wilcoxon rank-sum test). "a" and "b" indicate significance difference among sites (Mann-Whitney). T, water temperature; EC, Electrical Conductivity; TUR, Turbidity; TP, Total Phosphorus; TN, Total Nitrogen; DO, Dissolved Oxygen; BOD, Biological Oxygen Demand.
FIGURE 2 | Percent change of nutrients and BOD by natural riverine wetlands of southwest Ethiopia (Namely; Boye, Fisho, and Kitto wetlands) during dry and wet season. TP, total phosphorus; PO 4 , orthophosphate; TN, total nitrogen; DIN, dissolved inorganic nitrogen; BOD, biochemical oxygen demand (see Supplementary  Table S9 for measured parameters for each site during both the dry and wet seasons). Supplementary Table S3).

Composition and Diversity of Macroinvertebrates
A total of 57 macroinvertebrate families were identified during wet and dry seasons. Total taxa richness of Boye was significantly lower than Kitto (Mann-Whitney, p = 0.001). The minimum taxa richness of five was recorded in Boye stream (B2). The highest Shannon's diversity value was recorded in Kitto compared to Boye and Fisho wetland sites. The relative abundance of filter collectors (%FC), was significantly lower in Boye than Fisho, and Kitto wetland sites (Mann-Whitney, p < 0.01). No filter collectors were found in Boye wetland sites before joining the wetland (B2, B3, B4), except in B1 (Figure 3).
The Ethbios index was also lowest in Boye wetland sites. The incoming streams of Boye (B1, B2, B3, and B4) showed bad water quality and heavily degraded ecological status (Ethbios score 8-10). Most of the sites associated with Boye wetland depicted a major ecological disturbance (Ethbios score < 44). Fisho sites had a moderate water quality with significant ecological disturbance. The incoming streams of Fisho revealed moderate (Ethbios score = 58) and poor (Ethbios score = 20) water quality. Two sites (K1 and K5) associated with Kitto wetland, had a slight ecological disturbance (Ethbios score of 93 and 72, respectively) and some with significant ecological disturbance (Ethbios index score between 45 and 59) (see Supplementary  Table S4).
As shown in Figures 3A-C, the pooled value of sites before joining the wetland, within the wetland, and after passing through the wetland, presented a clear gradient in taxa richness, Ethbios FIGURE 3 | Diversity and compositional metrics of macroinvertebrates: (A) taxa richness, (B) Ethiobios index, (C) Shannon diversity and diatoms: (D) Shannon diversity-diatom, (E) taxa richness-diatom, (F) percent of tolerant species in sites before joining the wetlands ("before"), within the wetlands ("within"), and after passing through the wetlands ("after") in both dry (color: blue) and wet (color: red) seasons. The values are extracted from all the three wetlands i.e., "Boye," "Fisho," and "Kito" wetlands and the mean values are used.
Frontiers in Environmental Science | www.frontiersin.org index, and Shannon's diversity index values. The values of all three biotic indices were lowest in sites before joining the wetland.

Composition and Diversity of Diatoms
From both seasons, dry and wet, a total of 101 diatom taxa were identified. Diversity and composition varied among the three wetlands. A minimum of five and maximum of 29 diatom taxa was recorded in Boye (B5) and Kitto upstream (K1) respectively. In general, Boye sites showed significantly lower taxa richness than Kitto sites. Shannon's diversity index value was significantly lower in Boye than Fisho and Kitto wetland sites. Percent of tolerant taxa (% PT) was significantly different among the three wetland sites (Kruskal-Wallis, p = 0.001). Boye had the highest percent pollution tolerant taxa (p < 0.01) followed by Fisho (p < 0.01) and Kitto. A significant difference was observed between the three groups (Kruskal-Wallis, p = 0.001).
The pooled value of sites, before joining the wetland, within the wetland, and after passing through the wetland (Figures 3D-F), revealed a higher taxa richness of diatoms in sites after joining the wetland during both seasons. However, the taxa richness was higher during wet season than the dry season. The Shannon's diversity was the lowest and percent pollution tolerant taxa was the highest in sites before joining the wetlands than those within the wetland and after passing through the wetland.

Environmental Gradient and Macroinvertebrate Community Assemblage
In the dry season, the redundancy analysis (RDA) model selected NH 4 (p = 0.005), TN (p = 0.005), and DO (p = 0.035) as set of significant variables determining the distribution of macroinvertebrate families. The model explained 58.8% of the taxa-environment relation, where the first axis explained 29.7% and the second axis explained 24%. The overall model was significant with a p-value of 0.005 (Monte Carlo permutation test, n = 999). The first axis of RDA (Figure 4, Dry season) showed a pollution gradient. This axis was positively correlated with PO 4 and negatively correlated with DO. Streams feeding Boye wetland (B1-B4) were characterized by high PO 4 and TN and associated with the dominance of Culicidae (clu), Syrphidae (syr), Chironomidae (chi), and Oligochaetes (oli). In the negative first axis, sites of Boye located within the wetland (B5-B7) and at the outlet point (B8) were characterized by high turbidity. The dominant macroinvertebrates in these sites were Dytiscidae (dyt), Belostomatidae (bel), Corixidae (cor), Coenagrionidae (coa), Veliidae (vel), Hydrophilidae (hyd), Hirudinea (hir), and Notonectidae (not). Axis 2 (24%) showed a pH and NH 4 gradient. It separated Fisho (F) sites, except F5, along with high level of NH 4 and lower pH. In the dry season, the dominant taxa in this sites were Sphaeriidae (sph), Naucoridae (nau), Caenidae (ca), and Libellulidae (lib). Kitto sites (K2, K3, K4, and K6) were grouped together with F5. The upstream of Kitto wetland (K1) was characterized by low level of PO 4 , and K2, K3, and K6 are positively correlated with PO 4 .
During the wet season, NO 3 and TN were selected as the set of important variables associated with the distribution of macroinvertebrates (Monte Carlo permutation, n = 999, p = 0.005). Axis 1 (38%) and axis 2 (26%) of RDA (Figure 4, Wet season) showed nutrient and pollution gradient respectively. Boye sites were characterized by high dissolved inorganic nitrogen (DIN) and BOD, whereas, Fisho sites were negatively correlated with DIN and BOD. The first axis grouped Boye sites into two, those sites before Joining the wetland (B1, B2, B3, and B4), being in one group and sites within the wetland (B5, B6, and B7) and at the outlet point (B8) being in the other group. The first group (B1-B4) are positively correlated with NO 3 , BOD, TN, and DO. Unlike the dry season, Simuliidae was abundant in Fisho sites (highly abundant in F4) during the wet season. Potamonidae dominated the upstream (F1) of Fisho wetland.

Environmental Gradient and Diatom Community Assemblage
In the redundancy analysis (RDA) model of the dry season that was built using a forward selection with a maximum permutation of 999, the set of important environmental variables determining the distribution of diatom taxa were pH (p = 0.02) and EC (p = 0.02). The overall model was significant in explaining the taxa-environment relationship (p = 0.01). The RDA model explained 71% of taxa-environment relation by the first two axes.
During the wet season, TP, NO 3 , and NH 4 were identified as the set of important environmental variables determining the distribution of diatom taxa among sites (Figure 5, Wet season). Axis 1 (46%) showed a nutrient gradient. The RDA triplot grouped some of Boye sites (B1, B2, B3, and B6) toward high NO 3 -N and TP.

DISCUSSION
Quality of the Water Before, and After Joining the Wetlands In Jimma, the natural riverine wetlands of the town receive streams which are considered as major carriers of the towns  waste. These wetlands exhibited a potential to temporarily retain nutrients (TN, NH 4 , TP, and PO 4 ) and organic pollutants (BOD). During the dry season, water quality improvement was observed by the significant reduction in concentration of TN and DIN, and a considerably lower NH 4 after passing through the wetland. Concerning phosphorus, a retention of TP and PO 4 was revealed from lowered concentrations in sites within the wetland and after passing through the wetland. These results are in line with a review done by Fisher and Acreman (2004) on the potential of natural wetlands to remove nutrients. Although there was a slight increase in the concentration of NO 3 in sites located within the wetlands, a significant reduction was observed in the effluent. The increased concentration within the wetland can be attributed to the conversion of ammonia to NO 3 or limited denitrification (Gómez et al., 2002). Unlike the dry season, no clear gradient was observed during the wet season.
The integrative response of biological indicators revealed the water purifying effect of the studied wetlands and their potential to support biodiversity. A considerable increase in taxa richness and Shannon's diversity of macroinvertebrates was observed in sites within the wetland and after passing through the wetland compared to those before joining the wetland. The relative abundance of filter collectors is among the recommended indices to assess ecological status of wetlands in Ethiopia (Mereta et al., 2012). The %FC and Ethbios index were also higher within the wetland. It can be associated with the reduced load of nutrients and organic pollutants after joining the wetlands. This agrees with the negative correlation of %FC with TP and PO 4 , and Ethbios with TN (see Supplementary Table S7). A decrease in %FC with increasing level of organic pollution was also reported by Gebrehiwot et al. (2017) while evaluating the feeding interaction along a pollution gradient. Lougheed et al. (2008) revealed a higher diatom taxa richness in sites with a relatively better quality than the most degraded sites. In this study, the higher diatom taxa richness and lower percent pollution tolerant taxa in sites after joining the wetland indicated water quality improvement. Shannon's diversity of diatoms also showed an increasing gradient within the wetlands and at the outlet. Additionally, negative correlation of diversity indices with DO and EC, and positive correlation of %PT with TP, TN, and DIN was observed (see Supplementary Table S8). Therefore, the observed positive gradient of the indices after joining the wetlands clearly indicates the wetlands' potential to purify water and to support biodiversity.
During the wet season, in areas of low flow rate, a higher water depth provides an increased attachment surface for diatoms by covering stone surfaces with water (Bojorge-García et al., 2014). Bojorge-García et al. (2014) reported a higher abundance of diatoms during the rainy season in tropical mountain streams of Mexico. Similarly, our finding indicated a higher abundance and richness of diatoms during the wet season. On the other hand, diversity indices of macroinvertebrates were lower in the wet season. Gebrehiwot et al. (2017) also indicated a similar finding. This might be due to habitat instability, resulted from increased runoff, which in turn drift off macroinvertebrates.

Purifying Efficiency Among Wetlands and Ecological Quality
The percent change of nutrients and organic pollutants differ among wetlands. This variation can be attributed to the difference in influent quality, and extent of catchment disturbance. Among wetlands, the concentration of nutrients and electrical conductivity were significantly higher in sites associated with Boye wetland than Fisho and Kitto wetland sites. This is because the streams associated with Boye wetland are the major carriers of Jimma town's liquid waste (Ambelu et al., 2013). The more polluted Boye wetland showed a higher TP reduction during the dry season. This can be due to the higher phosphorus input to Boye wetland from the streams flowing to it or anthropogenic pressure in the surrounding since its located in urban downstream. The wetlands nutrient retaining effect was not observed in wet season, rather there was more leaching of nutrients than reduction. This might be due to instability of sediments and habitats.
The biological indices reflect major ecological degradation in Boye wetland. High load of organic and inorganic pollutant makes the environment less favorable for sensitive taxa (Eppink et al., 2004). The highest percent of tolerant diatom taxa and the lowest Shannon's diversity, macroinvertebrate and diatom, was observed in Boye than the other two wetlands. This can be attributed to the heavily degraded incoming water to Boye. Compared to Boye wetland the other two wetlands, located in urban upstream, still have the capacity to sustain their ecological quality. This was indicated by the higher diversity and richness of the assessed bioindicators.
The Relation Among Environmental Gradient, Biota, and Sites In this study, the RDA model reflected the dominance of macroinvertebrate taxa which are tolerant of heavy pollution, in the incoming streams of Boye (B1-B4) and characterized by a high load of nutrients (PO 4 , TN, and NO 3 ) and organic matter (BOD). Beyene et al. (2009) reported the absence of macroinvertebrates in heavily polluted streams of Addis Ababa, Ethiopia. However, in this study, heavily polluted sites recorded a minimum of five. This discrepancy may be due to the difference in the type of pollution since the sites are located in a different land use type. The presence of moderately tolerant macroinvertebrate taxa in Boye sites after passing through the wetland indicated a relatively improved water quality. Most of Fisho and Kitto sites were dominated by moderately tolerant taxa. Highly tolerant taxa (Lakew and Moog, 2015) such as Oligochete were among the dominant during the wet season. During the dry season, the upstream of Kitto (K1) is characterized by a high level of DO and had a low level of organic pollution. This variation in macroinvertebrate communities among sites can be explained by the difference in habitat, and in the level of degradation. In general, the relation among macroinvertebrate community assemblage, environmental parameters, and sites, discriminate Boye wetland sites as the most degraded in both seasons. A previous study (Ambelu et al., 2013) also pointed out the presence of heavy degradation on this site.
Diatom taxa which are tolerant of high level of EC, e.g., Craticula ambigua (Alakananda et al., 2011) and heavy organic pollution and nutrients, e.g., Navicula gregaria and Gomphonema parvulum var. parvulum (Kelly et al., 1995;Gómez et al., 2008) were abundant in Boye sites along a high level of pH and EC. Gomphonema parvulum is among the taxa which are potential indicators of eutrophication (Bere and Mangadze, 2014). Less tolerant taxa, e.g., Gomphonema augur (Bauer et al., 2007) were found in Kitto upstream site (K1), which indicates good water quality. There was a homogenized distribution of highly tolerant, and moderately tolerant taxa during the wet season. This could result from the drifting effect of runoff since it is the major factor governing the seasonal variation in tropical wetlands ecosystem (Bellinger et al., 2006;Mitsch et al., 2010).
In indicating the water quality improvement and ecological quality, macroinvertebrates and diatoms, indicated a coherent response toward pollution gradient during the dry season. Similarly, a consistent response of these bioindicators was reported by Resende et al. (2010) in assessing water quality of the UI river, Portugal. Both bioindicators discriminated Boye site as the most ecologically degraded.

CONCLUSION AND RECOMMENDATION
In conclusion, this study proved the potential of natural riverine wetlands in Jimma, Ethiopia, to retain nutrients and organic pollutants. The gradient of nutrients, BOD, and DO before and after joining the wetlands showed retention of pollutants within the wetlands. Percent change of nutrients and organic pollutants varied among wetlands and between wet and dry seasons. During the wet season, a higher concentration of nutrients was observed within the wetland. However, the water quality was better (diluted) during the wet season. Macroinvertebrates and diatom metrics clearly indicated the water purifying effect of studied wetlands and their potential to support biodiversity. B1 and B3 of Boye were identified as severely degraded, whereas, K1 (Kitto upstream) was the only site with an intrinsic good ecological quality. In general, the integrated response of biota and environmental parameters discriminate Boye wetland sites as the most ecologically degraded. Hence, we conclude that, when the incoming water is heavily degraded, the retention of pollutants seriously affects the ecological quality of the wetlands themselves. Therefore, future integrated wetland management interventions should also target the incoming streams or rivers to reduce the nutrient and organic loading.
Therefore, this study suggested that when the influents are marginally polluted the natural riverine wetlands can serve as natural treatment systems without affecting their ecological quality. Simultaneous use of macroinvertebrates and diatoms along with environmental parameters is recommended to have a comprehensive information on water quality and ecological status of natural riverine wetlands. Future studies should include wetland macrophyte identification, and hydrological parameters, e.g., hydraulic retention time, to identify the major processes governing the nutrient removal.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.

AUTHOR CONTRIBUTIONS
AS contributed in data collection, analysis, and manuscript writing. AA, AB, and LT contributed in designing the study, data analysis, and manuscript writing. IS contributed in manuscript writing.

FUNDING
This study was financially supported by VLIR-UOS and VUB (BAS42).