Monitoring abiotic influences on aquatic communities provides a powerful perspective for understanding current and future assemblage structure dynamics. Specifically, stream fish assemblages are influenced by numerous abiotic features, including but not limited to water quality parameters, stream size/order, local and regional habitat features, and disturbance regime/magnitude (Harrell, 1978; Matthews et al., 1988; Taylor, 1997; Wang et al., 2003; Mollenhauer et al., 2019; Stearman et al., 2019). Generally, environmentally benign habitats, which present relatively consistent flow regimes, temperatures, dissolved oxygen concentrations, and nutrient levels, support stable fish assemblages spatiotemporally (Matthews et al., 1988; Helms et al., 2009). In contrast, systems subjected to harsh environmental conditions, including large fluctuations or unpredictability in the aforementioned features, may contribute to lowered fish assemblage stability across time (Ross et al., 1985; Helms et al., 2009). In addition to underlying, background, abiotic variation within freshwater systems, stochastic disturbance factors may further modify fish assemblage structural dynamics in unpredictable manners (e.g., Helms et al., 2009; Geheber and Piller, 2012; Piller and Geheber, 2015).

Generally, disturbances (applicable to a wide variety of ecological systems) include any environmental fluctuations and/or events that alter population structure and influence the availability of resources (Pickett and White, 1985). However, the type and magnitude of disturbance often determines level of influence on fish assemblage structural dynamics, particularly the recovery potential of an assemblage following a specific disturbance event. Past works have documented impacts of floods (e.g., Harrell, 1978; Meffe and Minckley, 1987; Matthews et al., 1988), droughts (e.g., Larimore et al., 1959; Marchetti and Moyle, 2001; Magoulick and Kobza, 2003), and storm events (e.g., Van Vrancken and O'Connell, 2010; Geheber and Piller, 2012) on stream assemblage structure. In general, such events may negatively impact less tolerant species while allowing generalist species to persist and possibly take over habitats previously occupied by less tolerant species (Vázquez and Simberloff, 2002; Wilson et al., 2008). In addition to naturally occurring disturbance events, anthropogenic disturbances may further modify and alter stream conditions and, in turn, influence fish assemblage structure. For example, flow alterations due to anthropogenic stream impoundments have been shown to have dramatic impacts on fish assemblages by allowing the replacement of intolerant species with those that are able to acclimate to novel harsh environments, i.e., broadcast spawners replace simple nesters, active swimmers replace benthic fishes, etc. (Carlisle et al., 2011). Such human-caused disturbances may intensify effects of natural disturbances by either prolonging their duration or creating situations where a natural disturbance becomes amplified (Mann and Emanuel, 2006; Banholzer et al., 2014). Ultimately, increasing human densities (i.e., urbanization) has led to decreases in aquatic indices of biotic integrity as well as shifts in species assemblage structure (Paul and Meyer, 2001; Walsh et al., 2007). Increases in urban land cover have resulted in sensitive taxa decreases and tolerant taxa increases in fishes and macroinvertebrates (Morley and Karr, 2002; Morgan and Cushman, 2005), and importantly, this widespread rapid urbanization has led to greater wastewater discharge into aquatic ecosystems (Li et al., 2012).

Wastewater effluents may alter freshwater environments directly by elevating temperature and flow via thermal discharge (Kinouchi et al., 2006). Specifically, Kinouchi et al. (2006) found that water temperatures of Ara River (Japan) tributaries were significantly elevated due to heated effluents produced by nearby upstream WWTPs. Additionally, WWTPs have been shown to increase pH and decrease dissolved oxygen concentrations due to sewage inputs rich in organic carbon, nitrogen, and phosphorous (Gücker et al., 2006; Figuerola et al., 2012). Moreover, these WWTP input effects have been shown to have dramatic impacts on fish assemblage structure (Guidetti et al., 2003; Azzurro et al., 2010; McCallum et al., 2019). McCallum et al. (2019) found that WWTPs may act as ecological traps by warming the water, attracting stress-tolerant species (sometimes including tolerant nonnatives) to outfalls. Effluents from WWTPs may provide stable temperatures close to outfall sites during summer months, while significantly increasing downstream temperatures during winter months, which might influence fish movements towards or away from these thermally altered environments (McCallum et al., 2019; Mehdi et al., 2021). Moving downstream from WWTP outfalls, fish assemblages may shift away from stress-tolerant species (McCallum et al., 2019). Benthic fishes have been found to drive changes in assemblage structure between effluent-impacted and non-impacted sites (Tetreault et al., 2013; Mehdi et al., 2021), suggesting the need for documenting benthic assemblage structure within impacted systems. For example, Mehdi et al. (2021) found Neogobius melanostomus (Round Goby) occurred in higher abundances near WWTP outfalls compared to Percina caprodes (Common Logperch), which was significantly more abundant nearly 3 km downstream of a WWTP outfall in Hamilton Harbor (Canada). In addition to impacts of wastewater effluent on specific species, significant declines in species abundance and species richness have also been reported in wastewater-impacted stream sites compared to unaffected sites (Tsai, 1968; Galib et al., 2018). Despite known impacts of WWTPs on stream biota, there are some instances where fish populations have been shown to recover following efforts to improve wastewater treatment practices by separating the effluent process and cooling waters, as well as utilizing cleaner technologies (Dauba et al., 1997).

Currently, we lack a complete understanding of the specific impacts wastewater effluents may have on temperate stream fish assemblages. Therefore, quantifying existing effects of WWTP effluents on freshwater fish assemblages may be crucial for establishing current and future conservation efforts within impacted systems. Working within the upper reaches of Niangua River drainage (Ozark Highlands, Missouri, USA), and focusing on stream riffle habitat, our objectives were four-fold. First, we set out to examine potential influences of WWTP effluent presence on fish species richness (i.e., α diversity). Second, we tested whether the presence of WWTP effluent had significant impacts on fish assemblage structure (i.e., β diversity). Third, we compared abundance and distributions of Cottus bairdii (Mottled Sculpin) in both the presence and absence of wastewater effluent. Specifically, our third objective was aimed at illuminating potential impacts of WWTP effluent on distributions of heat-intolerant and pollution-sensitive fishes such as C. bairdii (Mebane, 2001; Adams et al., 2015). Finally, we compared abiotic variables between WWTP-impacted and non-impacted riffle sites to determine whether habitat properties typically associated with wastewater effluent were present. For the first three objectives, we hypothesized that the presence of WWTP effluent would influence fish assemblage diversity measures overall, as well as abundances of the thermally sensitive C. bairdii. For the fourth objective, we hypothesized that WWTP effluent presence would be linked to elevated temperature and discharge rates.

MATERIALS AND METHODS

Study system.—The Marshfield WWTP is located in northwestern Marshfield, Missouri, USA, north of Interstate 44. Serving 7,562 people as of 2019, the facility is an EPA and Missouri Department of Natural Resources certified complete mix activated sludge wastewater facility, providing treatment of 1.5 million gallons of sewage per day and as much as 3.5 million gallons of wet weather inflow per day to its collection system (Missouri Department of Natural Resources, 2018). Effluent from the facility empties into the West Fork Niangua River, which ultimately flows into the main stem of the Niangua River (Fig. 1). The total WWTP effluent discharge is 0.066 cubic meters per second (m³s^–1), which is added to the West Fork via two outfalls (Missouri Department of Natural Resources, 2018).

Fig. 1

Map depicting the 26 sites sampled in the East (n = 10, gray circles) and West (n = 16, black circles) Forks of the Niangua River during the study. The three circles with X's represent the sites where temperature loggers were deployed, and the WWTP outfall is located immediately south of W1 (right). Map inset (upper left) depicts the location of Missouri in the U.S., and inset (lower left) shows the location of the upper Niangua River drainage in Missouri.

The East and West Forks of the Niangua River, both fourth order tributaries, converge at the Niangua River main stem ∼11 km north of Marshfield (Fig. 1). The Niangua Watershed is located within the Ozark Highlands Plateau physio-graphic region, which is dominated by Ordovician and Cambrian rock, mostly consisting of dolomite (Schulz, 2001). Due to the subsurface movement of water through dissolved channels in the dolomite, karst features like caves, sinkholes/losing streams, and springs are abundant in the Niangua Watershed (Schulz, 2001).

For our study, stream site selections were based on availability of riffle access via road crossings, general site proximity to the Marshfield WWTP, and generating an “even” spread of sites across each of the two stream gradients of interest. The WWTP outfall is in the upper reaches of the West Fork Niangua, and therefore sites within the West Fork were considered effluent impacted (Fig. 1). All effluent-impacted sampling sites were located downstream of WWTP outfall due to lack of consistent water presence above the outfall. Sites along the East Fork Niangua were considered not impacted by effluent due to WWTP absence. Aside from WWTP presence in the West Fork, both forks were comparable in size (e.g., length, width, depth). Sixteen riffle sites were sampled in the West Fork, impacted, and ten riffle sites in the East Fork, not impacted (Fig. 1).

The Niangua River features a diverse assemblage of Ozark stream fishes, including numerous species of darters (Percidae) and minnows/shiners (Leuciscidae). However, an added focus for the present study were sculpin (Cottidae) due to their known sensitivity to environmental/heat degradation. The Niangua River system possibly harbors three congeneric species, including C. bairdii, C. carolinae (Banded Sculpin), and C. hypselurus (Ozark Sculpin). Although the latter two species are not well documented in the Niangua upper reaches where the present study was conducted. Generally, sculpins have limited home ranges, although they are known to disperse upstream, even in events of low flow (McCleave, 1964; Brown and Downhower, 1982; Greenberg and Holtzman, 1987; Wells et al., 2017). Given that the uppermost reaches of the West Fork Niangua River have been subjected to WWTP effluent input, assessing the abundances and distribution of sculpin may be important for understanding potential WWTP impacts on fish assemblages.

Fish assemblage sampling.—Between 25 May and 2 September 2019, fish sampling was conducted in riffle sites of the East and West Forks of the Niangua River. Each site was sampled once during this period (see Supplemental Table 1 for details; see Data Accessibility). All fish collections were conducted via kick seining using a weighted 1.8 m seine net (3 mm mesh) positioned downstream of the kicker, perpendicular to flow. In addition to summer sampling, we re-sampled four sites one additional time (two West Fork and two East Fork sites) during February 2020 (see Supplemental Table 1 for details; see Data Accessibility). Winter samples were not included in statistical testing due to lack of appropriate replication; however, these samples were included in multivariate visualizations to provide an anecdotal look at seasonal assemblage structure (i.e., summer vs. winter). For all sampling events, riffle sites were sectioned into 1.5 × 2 m kicking areas, and two kick seining passes were completed within each kicking area. All individuals captured were identified to species in the field, and individuals that could not immediately be identified in the field were euthanized in MS-222 (250–300 mg/L) and fixed in 10% formalin prior to identification and long-term storage in 70% ethanol. All preserved specimens were catalogued in the University of Central Missouri (UCM) ichthyology research collections.

Fish assemblage data analyses.—To compare α diversity between the East and West Forks, rarefaction (interpolation) and extrapolation (prediction) of species richness was quantified using the R package iNEXT (Hsieh et al., 2016). This approach allowed us to examine accumulation of species richness across sampling events in each fork and provided richness predictions for hypothetical additional sampling efforts. To examine assemblage structure among sites (β diversity), summer fish sample data were square root transformed prior to multivariate analyses to lessen any potential influences of overly dominant species and/or potential differences in sampling effort across all sites. A Bray-Curtis similarity matrix was used to determine assemblage relationships among all 26 sites (Bray and Curtis, 1957; Clarke, 1993), and non-metric multidimensional scaling (nMDS) was used to generate visual representation of the structural relationships in two-dimensional space. This procedure was completed using the metaMDS function of the vegan package in R (RStudio version 1.2.1335; Oksanen et al., 2019). It should be noted that an additional nMDS including the four winter site samples was run and is presented in results alongside the summer sample nMDS; however, these additional data were not included in any other analyses presented. To test for significant assemblage differences between the East and West Fork sites, a one-way permutational multivariate analysis of variance (PERMANOVA) was run using the Bray-Curtis matrix based on Type I sum of squares and 9,999 permutations of residuals under a reduced model. Additionally, site upstream distance (km upstream from the confluence) was included as a covariate in the model to account for possible influences of spatial auto-correlation. The PERMANOVA model was run using PERMANOVA+ in Primer 7 (Anderson et al., 2008). Similarity percentage (SIMPER) analysis was used to determine which species were main drivers of assemblage structure dissimilarity between the two forks. This analysis was based on the initial square root transformed Bray-Curtis similarity matrix and was employed to examine the differing contributions of each species, in terms of mean abundance per site, between the two forks (vegan package in R; RStudio version 1.2.1335; Oksanen et al., 2019).

Additionally, analysis of covariance (ANCOVA) was run using the aov function (RStudio version 1.2.1335; Oksanen et al., 2019) to assess potential differences in abundance of C. bairdii between the two forks. Here, site upstream distance (km) from confluence was the continuous explanatory variable (covariate), abundance of C. bairdii was the continuous response variable, and location (i.e., East Fork or West Fork) was the categorical explanatory variable. For this analysis, a significant interaction term would indicate different trends (i.e., different slopes) for abundance of C. bairdii along the East Fork gradient compared to abundances along the West Fork gradient. Additionally, regression analyses examining relationships between abundance of C. bairdii and upstream distance were run for each fork independently to further characterize within fork trends (lm function; RStudio version 1.2.1335; Oksanen et al., 2019).

Habitat sampling and analyses.—To quantify stream habitat characteristics at each site, discharge, maximum and average depth, and substrate size were measured. Discharge of each riffle was taken by setting a measuring tape across the stream width, perpendicular to flow, at the midpoint of each riffle site. The measured width was then divided by ten (i.e., to measure at equal distances ten times across the width), and water depth (cm) and velocity (m/s) were measured at each of the ten points across the width. Each flow velocity measure was taken using an OTT MF pro velocity meter set 60% below surface depth (to measure mean velocity in water column), and a 20-second recording interval was used to determine mean velocity at each point. Discharge rate in cubic meters per second (m³s^–1) was then calculated following the procedure of USGS (1982) for each riffle site. For eight West Fork sites and five East Fork sites, discharge could not be measured due to equipment failure. Additionally, the uppermost and lowermost width (m) as well as the length (m) of each riffle were collected for sampling area comparisons between the two streams. For each riffle, the two widths were averaged, and the resultant value was multiplied by riffle length to obtain general riffle surface area (these data were collected for 13 West Fork sites and nine East Fork sites). An independent t-test comparing riffle surface areas of the West and East Fork sites was performed to verify the assumption that site ‘sizes’ were comparable (i.e., not significantly different) between the two forks. As assumed, there was no significant difference in mean surface areas between the East and West Fork sites (df = 16, t = 1.08, P = 0.298; see Supplemental Table 2 for details; see Data Accessibility).

Substrate size of each riffle was collected using a modified version of the zig-zag procedure developed by Bevenger and King (1995). Substrate size was taken by walking two diagonal trajectories across each riffle (starting from the downstream edge) and measuring ten particles total (one particle was measured with each of five steps for the first diagonal trajectory [n = 5], and one particle was measured with each step during the second trajectory [n = 5]). The intermediate side of each substrate particle was measured using a substrate sizer (AL-SCI Field Sieve, West Trenton, NJ). Substrate was later categorized as gravel (2–16 mm), pebble (16–64 mm), cobble (64–256 mm), boulder (>256 mm), or bedrock (Compton and Taylor, 2013). Habitat variables were measured at the time of fish sampling for each site. Additionally, Hobo Water Temp Pro v2 Onset (U22-001) temperature loggers were deployed at two sites in the West Fork (one just below the WWTP outfall and one ∼2.5 kilometers downstream from the outfall) and one site in the East Fork (Fig. 1). Loggers recorded temperature from 15 January 2020 to 15 February 2020 at one-hour intervals, and loggers were deployed during winter months to examine potential ‘warming ability’ of WWTP effluent release.

To examine summer abiotic variation among study sites, principal component analysis (PCA) was run based on the aforementioned habitat data collected at each site. Principal component analysis included the 13 sites where all abiotic parameters were measured. Due to some missing data among sites (e.g., discharge was not collected from 13 of the sites), we deemed it appropriate to also run a second PCA which only included substrate composition data, as these data likely vary the least amount on a daily and seasonal basis. All 26 sites were included in the second PCA. Prior to each PCA run, data were normalized to make all variable units comparable. Principal component abiotic variable loading values were reported for both runs, and bi-plots (PC1 and PC2 depicted) were used to visualize abiotic associations among East and West Fork sites. PCAs were run using the built-in R function princomp (RStudio version 1.2.1335; Oksanen et al., 2019).

Long-term temperatures recorded from the two sites in the West Fork and one in the East Fork were plotted across time and visually compared. Long-term temperature is presented from 15 January to 15 February 2020. It should be noted that temperature data were also collected during summer 2019; however, due to temperature loggers being lost from certain field locations, summer data were limited and are not presented here.

RESULTS

Fish species richness comparisons.—Cumulative raw fish species richness in the East and West Forks was n = 9 and n = 13, respectively. The rarefaction species accumulation procedure demonstrated a slightly lower species accumulation curve in East Fork α diversity compared to that of the West Fork; however, 95% confidence intervals for richness extrapolation within each fork show some overlap (Fig. 2).

Fig. 2

Rarefied species richness within both the East and West Forks. Solid lines represent interpolation of species accumulation and dashed lines represent extrapolation of species richness if additional individuals were to be captured (gray = East Fork; black = West Fork). The black triangle represents observed species richness in the West Fork, and the gray circle represents observed richness in the East Fork. 95% confidence intervals are depicted for both accumulation curves.

Fish assemblage structure.—PERMANOVA indicated that fish assemblages of the East and West Forks were significantly different (P < 0.001), and this was determined after accounting for the significant effect of the upstream site distance covariate in the model (Table 1). nMDS ordination (2D stress = 0.174) illustrated a clear visual separation between sites of the two forks in terms of fish assemblage structure (Fig. 3, upper). Assemblages at sites resampled during winter showed the same relationships that were found among our summer samples (i.e., in ordination space; 2D stress = 0.177); East Fork winter samples did not differ from East Fork summer samples, and West Fork winter samples did not differ from West Fork summer samples (Fig. 3, lower). Overall, SIMPER analysis (summer assemblage data only) reported 58.74% average dissimilarity between East and West Fork assemblages. Etheostoma spectabile (Orangethroat Darter), E. flabellare (Fantail Darter), and C. bairdii showed the three highest average dissimilarity to standard deviation ratios (1.53, 1.73, and 1.27, respectively) between the two forks (Table 2). For reference, a species with a high average dissimilarity to standard deviation ratio (e.g., >1) is a likely strong contributor to overall assemblage differences. Specifically, E. spectabile occurred at greater abundances in West Fork sites, while E. flabellare and C. bairdii occurred at greater abundances in East Fork sites. In fact, E. flabellare was absent from all West Fork sites. Other species of notable interest were Noturus exilis (Slender Madtom), which was found in higher abundances in the East Fork, and Luxilus zonatus (Bleeding Shiner), Chrosomus erythrogaster (Southern Redbelly Dace), and Campostoma spp. (Stone-roller), which all occurred in greater abundances in West Fork sites (Table 2).

Table 1

PERMANOVA table of results indicating a significant difference in fish assemblage structure between East and West Fork sites (Pseudo-F = 10.481, df = 1, P = 0.0001) after accounting for the effect of the site distance covariate, which does influence fish assemblage structure (Pseudo-F = 7.813, df = 1, P = 0.00015). Stream distance data (between sites) were square root transformed prior to analysis to ensure that distance variation of East and West Fork sites were more comparable.

Fig. 3

nMDS based on Bray-Curtis similarity, depicting fish assemblage relationships among all study sites (West Fork: black triangles represent individual sites and the dark shaded hull encompasses sites; East Fork: gray circles represent individual sites and the light shaded hull encompasses sites). Upper panel depicts relationships among all 26 summer samples; lower panel depicts relationships among all summer samples and the four sites revisited during the following winter. Sites revisited are labeled ending with “w.”

Table 2

Similarity percentage (SIMPER) showing species contributions (Contrib%) to dissimilarity between East and West Fork sites. West and East Fork average abundance (Av.Abund) columns list the transformed mean abundances for species in both forks, ‘Diss/SD’ represents the ratio between the average dissimilarity and standard deviation, and the cumulative % (Cum.%) column lists the cumulative dissimilarity across species with a 90% cutoff. The average dissimilarity between the East and West Forks was 58.74%.

Prior to ANCOVA, the assumptions of covariate and treatment independence, and homogeneity of variance were tested and met. We then ran two separate ANCOVA models, one with and one without the interaction term included, and we compared the two models via ANOVA to determine model fit. This procedure showed that the interaction was not significant, and removing the interaction did not impact model fit. The final ANCOVA model determined that the relationship between site upstream distance (i.e., the covariate) and abundance of C. bairdii was significant overall (F = 6.494, df = 1, P = 0.018), but abundances between site locations (i.e., East or West Fork) did not differ (F = 0.052, df = 1, P = 0.822; Fig. 4). Since upstream distance explained abundance of C. bairdii overall, we decided to examine relationships between upstream distance and abundance of C. bairdii within each fork independently using linear regressions. Within the East Fork, abundance of C. bairdii decreased significantly with distance from the confluence, i.e., upstream distance (linear regression, F = 6.192, df = 1,8, P = 0.038, R² = 0.366; Fig. 4). Within the West Fork, abundance of C. bairdii showed no significant relationship with upstream distance from confluence (linear regression, F = 1.489, df = 1,14, P = 0.243, R² = 0.0316).

Fig. 4

Relationships between abundance of C. bairdii (per site) and upstream site distances (km) from the lower confluence (with the Niangua River proper). West Fork sites are depicted with black triangles, and East Fork sites are depicted as gray circles. Linear trend lines show relationships in each fork and further illustrate findings from ANCOVA.

Habitat structure.—In general, the PCA which included the 13 sites where all abiotic variables were collected indicated that maximum depth, average depth, and percent gravel were generally greater in West Fork sites (Fig. 5, upper). This finding was based on variable loading values associated with PC1 variation explained (42.9%, Table 3). Variation explained by PC2 (27.7%) was most related to discharge among sites, and moreover, West Fork sites generally showed greater discharge (Fig. 5, upper; Table 3). For PC variable loadings, values approaching |1| indicate strong association between the given variable and the particular PC in question.

Fig. 5

Principal component analyses illustrating abiotic relationships among sites in both the East and West Forks. Upper panel depicts the abbreviated dataset (i.e., only includes the 13 sites where all abiotic variables were collected). Measures of depth and % gravel best explained variation represented by PC1 (i.e., values increase towards the positive), and discharge and % pebble explained variation represented by PC2. Lower panel depicts substrate relationships among all 26 sites. PC1 explains variation in % gravel and % cobble among sites, whereas PC2 explains variation in % pebble and further % gravel variation among sites. For both panels, West Fork sites are represented by black triangles and are encompassed by the dark shaded hulls and East Fork sites are represented by gray circles and are encompassed by the lighter shaded hulls.

Table 3

Principal component loading values for seven habitat variables across 13 sites.

The PCA based on substrate data alone resulted in a PC1 that accounted for 39.4% of substrate composition variance. Specifically, PC1 showed a positive association with percent gravel and a negative association with percent cobble (Fig. 5, lower; Table 4). Generally speaking, percent gravel was greater in West Fork sites, while percent cobble was greater in East Fork sites (Fig. 5, lower). Principal component 2 (33.6% variance explained) was positively associated with percent gravel and negatively associated with percent pebble (Fig. 5, lower; Table 4).

Table 4

Principal component loading values for the four substrate categories measured across all 26 sites.

Temperature logger data, which provided a more accurate and consistent long-term picture of temperature in both forks, showed that temperatures in the West Fork were consistently higher than those of the East Fork between mid-January and mid-February (Fig. 6). During this time frame, the West Fork loggers recorded average water temperatures of 8.05°C directly below WWTP outfall and 7.95°C further downstream of the outfall. Overall, winter temperatures in the West Fork ranged between 4.0–12.0°C with an average of 8.1°C (i.e., based on cumulative data from both West Fork loggers). Winter temperatures in the East Fork ranged from 2.8–10.7°C with an average of 6.6°C (Fig. 6).

Fig. 6

Stream temperatures from one East Fork site and two West Fork sites measured 15 January 2020 through 15 February 2020. East Fork site temperature is represented by a solid gray line, and temperatures from West Fork sites are represented by a black dotted line for the WWTP “outfall” and a black dashed line for the further downstream site.

DISCUSSION

Past studies have demonstrated how wastewater effluents (e.g., industrial, agricultural, or sewage treatment) directly alter stream conditions such as temperature, stream flow, dissolved oxygen, and nutrient levels (Hamdhani et al., 2020). Short-term, these altered conditions may immediately allow effluent-tolerant fishes to inhabit territories/microhabitats within the impacted area, while effluent-intolerant species are displaced (e.g., Porter and Janz, 2003; Piller and Geheber, 2015). However, consistent and continued effluent inputs in a system (e.g., sewage treatment wastewater) may cause wide-reaching environmental degradation, leading to the extirpation of less tolerant species from entire stream reaches (Gafny et al., 2000). Even though it is known that aquatic species extirpations and/or displacements are potential repercussions of wastewater effluents, testing the “real-time” influences of effluents on aquatic assemblages is often problematic because pre-effluent assemblage data are not always available or do not exist. Since pre-WWTP assemblage data were not available for the West Fork in our study, we relied on the neighboring East Fork as a reference stream (i.e., a proxy for pre-effluent). Because of the shared species pool between the two forks, and the general similarity in physical size and broad features of the two forks, our approach for investigating WWTP effluent impacts on fish assemblage structure was quite reasonable. Using this approach, we observed significant differences between East and West Fork fish assemblages, and moreover, these results suggest that observed assemblage distinctions may be influenced by wastewater effluent inputs (e.g., leading to elevated stream temperatures) in the West Fork.

Assemblage diversity and structure.—Our finding of greater cumulative raw species richness in the West Fork compared to the East Fork was not a surprise considering more sites were sampled in the West Fork. However, rarefaction and extrapolation suggests that α diversity is likely quite similar within both forks based on observed 95% confidence interval overlap. Confidence interval overlap aside, it is worth noting that interpolated and extrapolated richness values were greater in the West Fork, and at the very least this suggests WWTP effluent has not led to drastic reductions in species richness among riffles of the West Fork; however, this result provides no insight into how diversity is structured among sites.

Analyses of fish assemblage structure revealed a clear and significant distinction between East and West Fork sites driven by several small-bodied benthic fishes. Considering small-bodied benthic fishes often lack the mobility to avoid localized disturbances due to their site attachment (e.g., Guidetti et al., 2003; Tipton et al., 2004), we suspect the benthic assemblage differences we identified are long-standing distinctions consistent with the presence of WWTP effluents. Our winter site revisits, although not robust in sample number, show repeated patterns of the assemblage structure observed in the summer, suggesting that riffle assemblage structure distinctions may be consistent across seasons as well. Tetreault et al. (2013) found small-bodied benthic fishes, such as Etheostoma caeruleum (Rainbow Darter), to display low mobility and experience consistent reductions in abundance downstream of wastewater effluent outfalls, while abundances of mobile, larger-bodied fishes (e.g., suckers and sunfishes) experienced abundance increases. Darter sensitivity to anthropogenic wastewater may at least partially explain the distributions of key benthic species in our study. One of the most striking disparities between the East and West Fork fish assemblages was the absence of E. flabellare in the West Fork, while it was relatively abundant in East Fork sites. The complete absence of E. flabellare in the West Fork suggests that environmental conditions, potentially effluent exposure, may have been too harsh for this darter species. Scarce literature exists on the specific effects that wastewater effluents have on population densities of E. flabellare, although polluted waters and siltation are known threats to this species (Gelder and Cribb, 2015). Interestingly, we did not find the same distributional pattern for all darter species collected during the study. Etheostoma spectabile was collected in both forks but consistently in greater abundances in West Fork sites (i.e., effluent present). Etheostoma spectabile is extremely successful at adapting to a wide range of environmental conditions (Blair, 1959), and our results suggest it may have a greater tolerance threshold for wastewater effluent impacts compared to other co-occurring benthic fishes.

Our third study objective focused on potential impacts of WWTP effluent on distributions of the known heat-intolerant and pollution-sensitive C. bairdii. Although we set out to examine distributions of C. bairdii in the presence and absence of effluent a priori, we coincidently found this species also to be a large contributor to overall assemblage distinctions between the two forks (i.e., based on SIMPER). Considering we predicted that C. bairdii would be negatively impacted by wastewater effluent presence due to known sensitivity, identifying the species as an important contributor to assemblage structure was not surprising. Specifically, when assessing the relationship between upstream site distance and abundance of C. bairdii (per site), we found a significant negative correlation in the East Fork but no such significant correlation in the West Fork. Additionally, the absence of C. bairdii in upper sites of the West Fork is a notable distinction between the two forks (i.e., C. bairdii was collected from every East Fork site) and may suggest a localized impact of WWTP effluent on sculpin directly downstream of the outfall. Meanwhile, the overall downstream increase in abundance observed in the East Fork may demonstrate the expected abundance-distance relationship in the absence of wastewater effluent. This is important to consider, as studies have shown that fish abundances tend to increase downstream under normal environmental conditions as habitat complexity increases (Gorman and Karr, 1978).

As stated previously, differences in abundances of C. bairdii and E. spectabile between forks were, at least partially, responsible for overall fish assemblage structural differences. Interestingly, we observed opposing trends in the abundances of these two benthic species between the two forks. On average, E. spectabile were nearly three times as abundant in the West Fork as in the East Fork (mean abundance in the West = 20.44 individuals [per site]; mean abundance in the East = 7.4 individuals [per site]). In contrast, average abundance of C. bairdii was greater in the East Fork than the West (mean abundance in the East = 3.6 individuals [per site]); mean abundance in the West = 2.56 individuals [per site]). It is possible that this opposing pattern in species abundances is linked to competitive displacement of one benthic species by (or with) the other. In the absence of environmental disturbance, C. bairdii may outcompete E. spectabile, as it has been documented that E. spectabile reduces habitat use in the presence of C. carolinae, a closely related species of Cottus (Taylor, 1996). Furthermore, species less tolerant of environmental alterations, such as C. bairdii, which is known to have narrow temperature requirements, may be subjected to increased interspecific competitive pressures in the wake of abiotic environmental alterations (Greenberg, 1988; Van Zuiden et al., 2016). Therefore, in the West Fork, wastewater effluent may have shifted the abiotic environmental conditions outside (or to the edge) of C. bairdii tolerance, and competitive pressures of E. spectabile may have followed. This scenario potentially explains why C. bairdii was absent from sites closest to the effluent release, while E. spectabile was relatively abundant. All said, the disparity in abundance of C. bairdii between the two forks was much less than what was observed for E. spectabile, so it is not likely that competitive pressures alone are responsible for the contrasting abundance patterns.

Temperature and abiotic site structure.—An initial assumption of our study was that wastewater effluent input would influence stream temperatures. Although we were not able to collect continuous stream temperature data across the two forks for the entire study duration, we were able to collect continuous stream temperatures during winter. Since we were interested in directly observing whether WWTP effluents increase stream temperatures, we took advantage of the perceived disparity between natural winter stream temperature and effluent input temperature to best capture any effluent-driven temperature increases. Overall, we uncovered a striking temperature distinction between the two forks; winter temperatures in the East Fork were consistently lower than the corresponding temperatures (based on time of measure) collected from the West Fork. Previous studies have shown that WWTP effluents cause perpetually elevated stream temperatures, and such increases were especially noticeable during winter and early spring months (Kinouchi, 2007). In the present study, temperatures collected just downstream of the WWTP effluent outfall were between 1° and 2°C greater than temperatures in the East Fork reference site during January and February 2020. Past studies aimed at quantifying the impacts of WWTP effluent release have also found success documenting elevated stream temperatures in relation to effluent outfalls during winter months (Brown et al., 2011; Mehdi et al., 2021), and we suspect the temperature differences observed between forks in our study system is at least partially the result of WWTP effluent release. That said, the upper Niangua River has numerous spring water inputs (i.e., cold water inputs), and it is feasible that such inputs could further influence the temperature differences observed between forks. Additionally, it is worth noting that the maximum weekly average temperature tolerance of C. bairdii is estimated to be 24.3°C (Eaton and Scheller, 1996). Although winter temperatures presented were well below this mark, average temperatures were 24.05°C in the West Fork and 23.89°C in the East Fork during 8–10 July 2019 (i.e., this was the only summer temperature data able to be recovered). Although anecdotal, these data suggest that these systems come close to reaching C. bairdii thermal maximum tolerance during summer months, and even a moderately small temperature increase may have the potential to impact fish occurrence.

Overall, we were interested in understanding impacts of WWTP effluent on fish assemblages; however, we also quantified abiotic site relationships in an attempt to understand underlying ‘natural’ distinctions in site-specific habitat features. Because the impact of WWTP presence on riffle substrate composition is likely minimal within our study system, we specifically aimed at quantifying site substrate composition to characterize natural abiotic relationships among sites of the two forks. Substrate size is a major driver of benthic fish presence, as different species use specific substrates for foraging, reproduction, and protection (e.g., Geheber and Frenette, 2016). Gravel substrate appeared in greater proportions in West Fork sites than East Fork sites, while cobble cover was found in greater proportions in the East Fork. Although no sites in our study were characterized by silt or sand substrates, it is worth mentioning that fish assemblage biomass is typically highest in areas where substrates are primarily composed of larger particles, i.e., gravel–cobble, as opposed to majority small-sized particles, i.e., silt-sand (Schlosser, 1982). Past studies document species-substrate associations among benthic stream fish assemblages, with different species utilizing different substrate sizes, which may reduce interspecific competition (Schlosser and Toth, 1984; Geheber and Frenette, 2016). Because species–substrate associations likely exist in the Niangua River system, it is plausible that the significant fish assemblage difference (between East and West Fork sites) is partially driven by substrate composition distinctions between the two forks. However, the considerable overlap in substrate composition among East and West Fork sites ultimately suggests riffle assemblage differences are minimally related to substrate composition differences. It is worth noting that reaches directly below the WWTP outfall in the West Fork did have high bedrock presence (based on KAH personal observation), though these sites were not sampled for fish or habitat due to accessibility issues. Typically, systems exhibiting higher discharge rates and/or events will have larger substrate since smaller particles get pushed downstream (Herbst and Cooper, 2010). Therefore, the overall larger substrate and lower baseline discharge observed in East Fork sites may suggest it is a flashier system, in terms of peak flow events, than the West Fork.

Conclusions.—We demonstrate significant riffle assemblage differences between sites in the presence of WWTP effluent and sites where WWTP effluent was absent. Assemblage differences were characterized by several benthic species distinctions, including consistent presence of E. flabellare and greater abundances of C. bairdii, a heat-intolerant species, in sites of the East Fork (i.e., wastewater effluent absent). Overall, findings suggest that fish assemblage structural differences between the East and West Fork sites could be linked to wastewater effluent release in the West Fork (i.e., which has increased water temperatures). Additionally, underlying differences in East and West Fork substrate composition may also play a role in the observed fish assemblage differences. Future studies on how wastewater effluents may impact benthic fish behavior, physiology, and life histories may be beneficial to a better understanding of how these spatially limited populations respond to environmental disturbances (Hill and Grossman, 1987). It is also important to consider how natural disturbances may exacerbate WWTP effluent impacts, such as the effects of reduced water flow due to drought which may reduce dilution of pollutants attributed to sewage (Gasith and Resh, 1999). Through the present study we provide a greater understanding of how wastewater, perhaps combined with natural variation, may affect fish assemblages that are already spatially limited and therefore especially vulnerable to temperature elevations and pollutants associated with such effluents.

DATA ACCESSIBILITY

Supplemental material is available at  https://www.ichthyologyandherpetology.org/i2022089. Unless an alternative copyright or statement noting that a figure is reprinted from a previous source is noted in a figure caption, the published images and illustrations in this article are licensed by the American Society of Ichthyologists and Herpetologists for use if the use includes a citation to the original source (American Society of Ichthyologists and Herpetologists, the DOI of the Ichthyology & Herpetology article, and any individual image credits listed in the figure caption) in accordance with the Creative Commons Attribution CC BY License.