Stream baseline conditions shape functional responses to wastewater: evidence from insect-dominated sites

Macroinvertebrate data and study design

Monitoring data used in this study were derived from the Hessian Agency for Nature Conservation, Environment, and Geology, covering the Rhine-Main metropolitan area in the south and rural areas in the north of Hessen, Germany (see Enns et al., 2023 for details). The larger river catchments covered in this study are the Rhine, the Main, as well as the Weser and their respective tributaries (map of study region in Enns et al. (2023); WWTP overview in Table S3). Monitoring data encompassed macroinvertebrate assemblage data, which were collected using standardized WFD methods (Haase et al., 2004) from 2006 to 2014, with additional samplings conducted in 2017 and 2018. The sampling season for all sites ranged from February to July. The taxonomy of the monitoring data follows the German operational taxon list (Haase, Sundermann & Schindehütte, 2006), which defines the minimum identification requirements for Germany under the WFD and includes classifications at varying taxonomic levels, primarily from family to species.

To assign monitoring sites to each WWTP, an upstream and downstream site for each WWTP was selected according to the criteria described in Enns et al. (2023). In total, 169 site pairs located above and below WWTPs were selected for our main dataset. It is important to note that 34 sites had an intermediate location between a possible upstream and a possible downstream site. These intermediate sampling sites were treated both as an upstream site of one pair and a downstream site of the other pair. This was necessary in order to not further limit the dataset. In cases where site pairs were sampled in multiple years only the most recent year was selected for the analysis. Within site pairs, the average time between sampling dates was 6.5 (± 19.3 SD) days. The minimum distance between a WWTP and the respective downstream sampling site was 70 m and the maximum distance was 27,356 m. WWTPs are characterised in Table 1, showing that most are grouped into size classes two and four, with average population equivalents of 2,765 and 34,754 human inhabitants respectively. Most of the analysed WWTPs employ a mechanical and biological treatment with nitrification, denitrification, and phosphate precipitation (Fig. S1). Regularly monitored effluent parameters were available and included biological oxygen demand, total ammonia, and total phosphorus concentrations, which tend to increase with decreasing size class, indicating that smaller WWTPs tend to emit higher organic pollutant loads in our dataset.

Trait selection and extraction

In order to more accurately assign trait data to each taxon in our dataset, taxa were validated beforehand using the taxa validation tool from freshwaterecology.info (Schmidt-Kloiber & Hering, 2015). The trait data, extracted from freshwaterecology.info, comprised mostly traits derived from Tachet et al. (2010). For our analyses, nine trait groups representing 49 trait modalities were selected. We selected traits that were likely to be impacted by WWTPs through colmation of the interstitial, changes in flow and temperature regimes, and through the introduction of emerging contaminants, including: aquatic life stage, dispersion, feeding group, locomotion and substrate relation, reproduction cycles, reproduction modes, resistant forms, respiration, and fecundity (Arenas-Sánchez et al., 2021; Descloux, Datry & Usseglio-Polatera, 2014; White et al., 2017). While the traits are not explicitly designed to detect chemical pollutants, they represent broader ecological responses to the combined impacts of wastewater discharges. In cases where taxa had multiple aquatic life stages (i.e., mostly members of the Coleoptera), trait data were always assigned at the adult level. We used a trait gap filling approach to ensure good trait coverage across our taxonomic dataset (see Haase et al., 2023), with trait aggregation having been found to be a good proxy for expert trait allocation (Kunz et al., 2022). This resulted in a trait coverage of 100% across all taxa, with 32.4% of traits coded at the original taxon level, 1% at the sub-species level, 44.5% at the genus level, and 22.2% at the family level. Taxa with insufficient trait coverage were raised to the next highest taxon, for which sufficient coverage was available, resulting in 444 final taxa. Trait data were fuzzy coded and affinity scores were transformed into relative frequencies. From the taxa-trait matrix, a Gower dissimilarity distance matrix was constructed and subsequently used to perform a principal coordinate analysis (PcoA) in the ade4 R package (Thioulouse et al., 2018). To determine the occurrence probability of different trait combinations, we used kernel density estimations (Chacon & Duong, 2018), calculated on two principal components of the functional space which explained 25.6% of variance.

Community weighted means

Community-weighted means (CWM) are used to quantify the average value of a set of traits within a community, weighted by the abundance of each species. In other words, CWMs weigh each trait modality, of the relative frequency trait matrix, by the abundance of taxa assigned to that trait and averages them for each community. We used the FD R package (Laliberté, Legendre & Shipley, 2014) for the calculation of the CWMs. Delta CWM values, calculated as the upstream CWM minus the downstream CWM for each site pair, were tested against zero using one-sample Wilcoxon tests.

Alpha diversity indices

We calculated commonly used trait-based metrics measuring functional alpha diversity for all sites as well as functional beta diversity between corresponding site pairs. Due to the constraint of the number of taxa on functional components (Villéger, Mason & Mouillot, 2008), functional richness (FRic) was calculated using three functional components (m), explaining 36,1% of variance, while functional evenness (FEve) and functional dispersion (FDis) were calculated using the entire trait space. The calculations were performed using the FD R package following the code of Múrria, Iturrarte & Gutiérrez-Cánovas (2020). The optimal number of functional components was estimated by calculating the mean squared-deviation between the initial functional distance and the standardised distance in the functional space. FRic was calculated as the hypervolume of a convex hull in functional space (Villéger, Mason & Mouillot, 2008). FEve was calculated using the minimum spanning tree method (Villéger, Mason & Mouillot, 2008) and is a measure of how regularly taxa are distributed within the functional space. FDis was calculated as the mean distance of every taxon to the centroid of the community in functional space (Laliberté & Legendre, 2010). We further calculated functional redundancy as the quotient of Rao’s quadratic entropy divided by Simpson diversity of a community (Ricotta et al., 2016), using the adiv R package (Pavoine, 2020). Functional redundancy of upstream sites was tested against downstream sites via a paired Wilcoxon test.

Beta diversity indices

Functional beta components were calculated using the mFD R package (Magneville et al., 2022), whereby a different functional space was built and beta diversity was calculated using three functional dimensions explaining 28.2% of the variance. Functional beta diversity metrics measure the dissimilarity of occupied functional space between two or more communities. It can therefore be further decomposed into functional turnover, which indicates changes in functional composition, and a nestedness-resultant component, measuring to what degree communities resemble functional subsets of each other (Villéger, Grenouillet & Brosse, 2013). A high functional beta diversity can therefore either result from high turnover, meaning that communities are differently structured in regards of their functional strategies, or from a high nestedness-related component, indicating that functional strategies hosted by one community represent a small subset of another community. By deducing the intensity of the different components, one can infer if niche differentiation (in case of high turnover), niche filtering (in case of high nestedness), or functional convergence (in case of low functional beta diversity) drives the observed functional patterns between two or more communities (Villéger, Grenouillet & Brosse, 2013).

Null models

To account for the intrinsic relationship between taxonomic richness and functional diversity metrics (functional diversity is constrained by taxonomic diversity; Gotelli & Graves, 1996), we calculated null models by shuffling species names within the trait matrix through 999 randomisations (Swenson, 2014). This way, we kept the number of taxa and their abundances within each community constant, only changing their functional categorization. Trait associations were shuffled using the total taxa pool so as not to limit the number of possible randomisations to a point where statistical deductions became impossible. Calculation and analysis of null models for beta functional diversity metrics followed Villéger, Grenouillet & Brosse (2013). To be able to compare metrics between communities, standardised effect sizes (SES) were calculated by subtracting means of the null distribution from the observed values and dividing this value by the standard deviation of the null distribution (Múrria, Iturrarte & Gutiérrez-Cánovas, 2020; Baker et al., 2023). To determine if the observed SES values for each site (alpha diversity metrics) and site pair (beta diversity metrics) were outside the null distribution (a test of statistical significance and an indicator of a deterministic effect; Larsen et al., 2024), a two tailed test was used. Subsequently, we calculated the proportion of sites and site pairs with significantly different observed SES values to the null distribution. Differences in alpha diversity metrics, SES values between upstream and downstream sites, and between components of beta diversity SES were tested using a Wilcoxon test. Further, beta diversity SES were correlated with number of connected households, population equivalents, as well as log transformed biological oxygen demand, ammonia, and total phosphorus concentrations using a Spearman rank correlation coefficient test. These effluent parameters are regularly monitored by the relevant authorities and were therefore the only reliable and accessible indicators.

EPT dominant site analysis

Sites of good ecological status are prime targets for conservation practitioners and EPT taxa are generally viewed as sensitive towards pollution (Sinclair et al., 2024). To determine whether sites with a higher dominance of EPT organisms were more vulnerable to the impacts of WWTPs, we generated a subset of sites, keeping only those site pairs that had a proportion of at least 50% EPT individuals in the upstream community (hereafter EPT dominated sites). This resulted in 60 communities from 30 site pairs for which the same analysis as in ‘Trait selection and extraction’ was performed. The functional space was reconstructed based on these communities and functional alpha metrics were calculated using 8 dimensions, explaining 69.9% of variance and beta metrics were calculated using three dimensions, explaining 30.4% of variance. A brief taxonomic characterisation of subset of sites is given in (Section 2.1S).

Functional trait changes of the communities

Functional changes differed depending on whether we considered the complete dataset or the EPT dominated dataset. In the EPT dominated dataset, we found clear differences in the trait characteristics of macroinvertebrate communities up- and downstream of WWTPs (Fig. 1A), with differences becoming less evident when we considered all sampling sites together (Fig. 1B). At EPT dominated site-pairs, we found significant changes in trait combinations related to reproduction (voltinism, fecundity, reproduction technique), dispersal (aquatic life stages, dispersion technique, locomotion), and feeding (feeding type), with all these traits significantly deviating from zero (Fig. 1A).

Traits with a positive shift, meaning an increase in the richness and/or abundance of species with these traits, were ovoviviparity, interstitial locomotion, shredding and absorbing feeding types, aquatic passive dispersal, aquatic adults (i.e., hololimnic life cycle), and low fecundity (<100). The observed shifts in trait diversity, both positive and negative, indicate which traits are being filtered by WWTPs and point toward communities composed of more generalist organisms and fewer EPT taxa. Organisms with the aforementioned trait combinations often gain competitive advantages under certain stress conditions, compared to those that are less adaptable or more sensitive (Paz et al., 2022). For example, ovoviviparous organisms carry their eggs and potentially regulate certain environmental conditions, e.g., avoidance behaviour or oxygen regulation through gill fanning. Further, interstitial locomotion might be more beneficial compared to epibenthic modes of locomotion like crawling or swimming, since the organism is able to hide under gravel during high flow events, e.g., storm overflows of WWTPs (Marino et al., 2024). In addition, organisms with absorption feeding habits might profit from the increased nutrient loads typical of wastewater effluents.

Traits with a negative shift between up- and downstream sites were univoltinism, aquatic egg stage, eggs as resistant forms, cemented and isolated eggs, high fecundity (between 1,000 and 3,000), aerial active dispersion, and scraping feeding types, with the latter five showing generally high affinities in Ephemeroptera and Plecoptera, but also some Trichoptera, Coleoptera and Gastropoda (Fig. S2). In particular, the Ephemeroptera and Plecoptera are sensitive toward trace metals and municipal wastewater (Let et al., 2022), and thus WWTPs may be influencing environmental filters which select against these traits. For example, insects with aerial active dispersion often emerge during the late evening or at night and suffer from light pollution (Perkin et al., 2014). It is thus reasonable to assume that WWTPs do not impact this trait directly, yet their spatial proximity to settlements and the resulting increased light pollution does. Further, nutrients introduced by effluents would promote the growth of epibenthic algae, thus favouring organisms with scraping feeding types. However, since this trait is associated with many ephemeropteran, trichopteran, and coleopteran taxa, the feeding type appears to be only passively selected against by WWTPs. Péru & Dolédec (2010) found that, in terms of biodiversity and functional metrics, river sections downstream of WWTPs deviate stronger from minimally disturbed reference conditions compared to sites located upstream of WWTPs. Similar to our results, traits with the strongest contribution to the deviation were related to reproduction types as well as locomotion and substrate relation, with the former potentially being attributed as a coping response against effluents and their impacts (e.g., organic pollutants, sedimentation, etc.). Furthermore, Mor et al. (2019) showed that traits such as small body sizes, semivoltine life cycles, eggs as aquatic stage, eggs as resistance forms, and gill respiration decreased downstream of WWTPs. This partially coincides with our results, with the only difference being that, in our results, the trait ‘gill respiration’ is not changing and semivoltine lifecycle is increasing, though this deviation may have resulted from differences in analytical methodologies. It could be expected that an increase in shredders and a decrease of scraping feeding types results in increased litter break down rates and a stronger biofilm formation. While Pereda et al. (2020) report these expected patterns, Burdon et al. (2023) and Englert et al. (2013) observed decreased invertebrate driven decomposition, with a comprehensive meta-analysis by Brauns et al. (2022) also showing lower leaf litter decomposition and a higher net ecosystem production in wastewater impacted streams. Nevertheless, the changes in CWMs observed in our study suggest that sites with a higher proportion of EPT taxa, and by extension their unique functional roles, are altered by downstream WWTPs, corroborating a previous study by Enns et al. (2023).

When considering the entire dataset of all 169 site pairs, trait compositional shifts downstream of WWTPs were much less apparent (Fig. 1B). Here, only traits related to fecundity and locomotion significantly deviated from zero. However, the absolute difference in average abundance between upstream and downstream sites was generally higher when we only considered the EPT dominant subset (absolute difference in average abundance for taxa with prevalence higher than 50% in EPT dominated subset: 16.8; in complete dataset: 3.7). This could explain the lower deviation of CWMs, since they are influenced more by highly abundant taxa (Péru & Dolédec, 2010). Another factor is the increased heterogeneity of communities of the complete dataset, which contributes to the indistinctness of the CWM response. This suggests that WWTPs tend to have a more substantial impact on the functional diversity of communities where insect taxa comprise a majority and further shows how results of an upstream-downstream comparison can differ depending on the baseline (upstream) community condition. In general, interpretations of changes in CWMs are not always straightforward, since stressors can select directly or indirectly for or against traits. This circumstance is detailed in Verberk, van Noordwijk & Hildrew (2013), whereby traits are connected to each other by the species that are associated with them. Thus, without a good understanding of species autecology and accompanying environmental and habitat data, it is difficult to disentangle the mechanisms behind how WWTPs filter traits and by extension the species to which they are associated.

Changes measured on alpha diversity scale

In both datasets, no functional alpha diversity indices showed a significant difference between sites up- and downstream of WWTPs (Fig. S4). Functional redundancy is often regarded as a valuable metric for conservation, since it gives information on the vulnerability of functional strategies (Mouillot et al., 2014; Ricotta et al., 2016). Thus, stable functional redundancy between site pairs can also be viewed positively. Contrarily, Péru & Dolédec (2010) as well as Mor et al. (2019) reported lower functional diversity below WWTPs, the latter showing that the impact magnitude depends on streambed substrate.

Changes measured on beta diversity scale

Most site pairs showed no significant deviations from null model expectations, suggesting niche convergence rather than differentiation. In other words, the set of functions fulfilled by communities up- and downstream of the WWTPs are, for most site pairs, quite similar. This finding was consistent when we considered the EPT dominant sites separately and when we repeated our analyses using complete data from all 169 site pairs. This therefore suggests that WWTPs are generally not changing the composition of the functional strategies of invertebrate communities we analysed. However, when considering all 169 site pairs, 14% of the site pairs (23 pairs) deviated significantly from the null model expectations, while only 10% of site pairs (three pairs) deviated significantly in the EPT-dominated site subset. Further, standardised effect sizes of functional beta components showed no significant correlation with WWTP characteristics from Table 1 (Fig. S3). Thus, it is reasonable to assume that in most cases (86%), the impact of WWTPs does not lead to detectable changes in the functional strategies of macroinvertebrate communities. However, this should not be interpreted as evidence that WWTPs are incapable of driving functional change. Rather, such effects may remain undetected due to current data limitations—particularly the lack of traits directly linked to chemical sensitivity—and the potential buffering capacity of functional redundancy. The low observed functional beta diversity hints at regional environmental pressures, as opposed to local factors like WWTPs, homogenising functional strategies across assemblages. Evidence for strong local factors influencing the functional composition of communities comes from intermittent rivers, where intermittency is a strong driver of functional turnover in invertebrate communities (Aspin et al., 2018; Piano et al., 2020; Viza et al., 2024). There is no evidence to date that chemical stressors—even when present in low concentrations—cause similarly strong effects on community functionality (Alric, Geffard & Chaumot, 2022). Thus, for intermittent rivers, WWTPs can even have a stabilizing effect on functionality since they guarantee a more regular water supply (Brooks, Riley & Taylor, 2006).

Functional diversity: limitations, interpretation, and future considerations

Functional diversity measures are a valuable tool to infer and understand important ecological circumstances. However, two major shortcomings exist thus far. First, if the environment selects for species with certain traits, other traits indirectly get selected for because they are intrinsically associated to that species (and not the other way around). Second, complex environmental pressures might select for a certain combination of traits, rather than a single trait. Thus, trait patterns are context dependent and result in low discriminatory power of diversity metrics, potentially hindering the interpretation of causalities between stressors and observed trait patterns (Heino, Schmera & Eros, 2013; Verberk, van Noordwijk & Hildrew, 2013). For our study area and period, access to comprehensive environmental data was limited—particularly with regard to measurements of chemical loadings in both the effluent and receiving waters. In addition, the other potentially confounding factors, such as changes in the surrounding land-use practices (Vandewalle et al., 2010) and the condition of the in-stream and riparian habitats (Calapez et al., 2021; Sargac et al., 2021), could not be assessed. Nevertheless, these problems could be mitigated by an a priori selection of relevant traits and/or by standardizing relevant environmental factors across study sites (Statzner & Beche, 2010). The former is a relatively easy task, achievable by reviewing existing literature. However, data on environmental factors are often either missing and/or spatiotemporally incongruent with the invertebrate sampling sites (Santos et al., 2021).

While functional diversity indices can give insights into the status of traits and functional strategies of communities, from which theoretical conclusions can be drawn, the relationship between indices and in-field functional processes remains unclear. There are several factors constraining the relationship between ecosystem functions and functional diversity. First, trait data are inherently generalized and often aggregated at higher taxonomic levels, which can obscure species-specific and intraspecific variability (Bolnick et al., 2011). Moreover, trait expression is influenced by complex, context-dependent interactions with environmental conditions, including adaptive responses that may alter trait expression in response to local stressors (Jourdan et al., 2024; Jourdan et al., 2019); (Hernández-Carrasco et al., 2025). Second, the relationship between functional diversity and ecosystem processes is constrained by the non-independence of traits and abiotic factors (Hamilton et al., 2020). Litter decomposition, as an example, is a highly relevant process in freshwater ecosystems and has a direct link to food and feeding traits of invertebrates. Studies suggest that organic matter breakdown in streams is either not correlated with invertebrate functional diversity indices and breakdown related traits (Vos & Schäfer, 2017), or it is negatively related with functional richness (Rideout et al., 2022), showing that we are far from a complete understanding of the relationships between ecosystem processes and functional diversity. Further in-field measurements of community composition and ecosystem processes, together with mesocosm experiments, would provide valuable data, which, through new analytical tools, could bridge the gap between theory and practice for many freshwater ecosystem functions.