8 Discussion and Recommendations

knitr::opts_chunk$set(fig.path = "fig/0600-discussion/", dev = "png")

The results presented above paint a consistent picture: ecological condition on the Neexdzii Kwa mainstem improves from downstream to upstream, with a clear biological signal of mild nutrient enrichment near Houston that fades with distance upstream. Below, we interpret these patterns in the context of the watershed’s land use and water quality history, compare with previous monitoring, and outline recommendations for ongoing stewardship.

8.1 Site Gradient and Nutrient Enrichment

The three Neexdzii Kwa mainstem sites exhibited a consistent downstream-to-upstream gradient in benthic community condition. Multiple independent lines of evidence — taxonomic composition, functional feeding group structure, tolerance-weighted biotic indices, and multivariate ordination — converged on the same pattern: BUL-01 showed biological signatures of mild organic enrichment, BUL-04 was intermediate with high variability, and BUL-05 supported a community characteristic of excellent water quality.


At BUL-01, downstream of Houston, the community was characterized by elevated Chironomidae (24%) and Oligochaeta (7%), high relative abundance of collector-filterers (31%, almost entirely Hydropsychidae), and a mean HBI of 4.14 (range 3.91–4.31 across replicates) — a “very good” rating with possible slight organic pollution. These patterns are consistent with increased nutrient loading and fine particulate organic matter availability — Hydropsychidae are net-spinning caddisflies that thrive in nutrient-enriched reaches where suspended food particles are abundant (Barbour et al. 1999). Notably, Hydropsyche was the single most abundant taxon at BUL-01 (20% of total abundance), a pattern not observed at the other sites. The co-dominance of Hydropsyche with collector-gatherers (chironomids and ephemerellids) and the relatively low representation of sensitive specialists (shredders comprised only 2% of individuals) together suggest a community shifted toward nutrient-tolerant generalists.


BUL-01 integrates the cumulative upstream pressures present in the watershed between BUL-05 and the river mouth: agricultural and cattle rangeland use through the valley, the Knockholt Landfill, and the Houston wastewater treatment plant outfall 500 m upstream. The Houston WWTP is the closest and most proximate point source and accounts for the incremental step-change between BUL-04 and BUL-01, but the benthic signal at BUL-01 reflects a watershed-scale integration of pressures rather than attribution to any single driver.


This biological signal aligns with the historical water quality record. Total phosphorus at mainstem stations consistently exceeds the CCME eutrophication trigger, with elevated concentrations at the Houston STP reach. Remington and Donas (2000) documented periphyton biomass near Houston averaging 145 mg/m² chlorophyll-a with increasing filamentous green algae — exactly the conditions that would support the filter-feeding Hydropsychidae assemblage observed at BUL-01.


At BUL-05, the upstream site, the community was dominated by Lepidostoma (42% of abundance), a shredder-herbivore caddisfly associated with intact riparian function and allochthonous leaf litter inputs. EPT (mayflies, stoneflies, and caddisflies) relative abundance (84%), taxonomic richness (45 taxa), and a mean HBI of 2.58 (range 2.40–2.72, “excellent”) all indicate high ecological condition. The high proportion of specialized feeders (shredders 42%, scrapers 13%) at BUL-05 contrasts sharply with the generalist-dominated community at BUL-01 — specialized feeders are considered among the most sensitive indicators of stream health (Barbour et al. 1999).

8.2 Diversity and Dominance

Shannon diversity was higher at BUL-01 (H’ = 4.32, log2) than at BUL-05 (H’ = 3.59), despite BUL-05 having greater taxonomic richness (45 vs 37 taxa). This apparent paradox reflects the strong numerical dominance of Lepidostoma at BUL-05, which depressed evenness. High diversity indices at moderately impaired sites are a well-recognized limitation of diversity metrics in bioassessment — a community with many tolerant taxa distributed evenly scores higher than a species-rich community dominated by one or two sensitive taxa (Rosenberg and Resh 1993). This underscores the importance of using composition-based metrics (% EPT, % Chironomidae) and tolerance-weighted indices (HBI) alongside diversity measures when assessing ecological condition.

8.3 Multivariate Community Separation

PERMANOVA confirmed that community composition differed significantly among sites (Table 7.1), with homogeneous multivariate dispersion (Table 7.2) confirming that this result reflects true compositional differences rather than unequal within-site variability. The NMDS ordination (Figure 7.1) visually confirmed this separation, with replicate samples clustering tightly within sites and clear separation among site centroids.


The pattern of within-site variability was itself informative. BUL-05 showed the lowest multivariate dispersion (mean distance to centroid = 0.14), reflecting a stable, uniform habitat where Lepidostoma dominated all three replicates similarly. BUL-01 was intermediate (0.21) — the nutrient enrichment signal (Hydropsychidae, Chironomidae) was present in all replicates but varied in relative intensity, suggesting a pervasive but non-uniform influence across the reach. BUL-04 showed the highest within-site variability (0.26), with individual replicates ranging from near-reference (high EPT, low Chironomidae) to moderately impaired composition. This pattern is characteristic of a transitional zone, where fine-scale differences in substrate, velocity, or organic matter availability among kick-net locations produce a mosaic of micro-conditions and correspondingly patchy communities. This contrasts with both ends of the gradient, where conditions are more uniform — consistently enriched at BUL-01, consistently reference-quality at BUL-05.

8.4 Comparison with Historical CABIN Data

The CABIN open data archive contains benthic records from site MOR37 (corresponding to BUL-01) from two sampling events: 2004 (BC MOE-FSP Skeena Region study) and 2018 (BC-Wet’suwet’en ESI study). Both historical samples showed Ephemeroptera-dominated communities with high EPT relative abundance. The 2025 BUL-01 community showed a notably different structure: lower EPT, substantially higher Trichoptera (dominated by Hydropsychidae), and elevated Chironomidae.


These differences should be interpreted cautiously. Although all three sampling events were within 49 m of each other (confirming the same reach), seasonal timing varied substantially — from August 13 (2018) to September 10 (2004) to October 3 (2025) — and differences in invertebrate phenology across this window could influence community composition. The two historical samples were collected under different CABIN studies with potentially different subsampling and taxonomic resolution protocols, and each is a single unreplicated sample — limiting the strength of temporal inferences.


Despite these caveats, the 2025 community composition at BUL-01 — with its high Hydropsychidae abundance and elevated Chironomidae — is consistent with a site experiencing mild nutrient enrichment. Whether the compositional differences between years reflect a real temporal change in site condition, interannual variability, or methodological differences can only be resolved through continued standardized monitoring at this location.

8.5 Implications for Restoration Planning

These results provide a spatially explicit benthic bioassessment of the Neexdzii Kwa mainstem and establish a baseline against which future changes can be measured. The clear site gradient — from mild enrichment at BUL-01 to excellent condition at BUL-05 — provides an empirical framework for prioritizing restoration and monitoring activities. Key implications include:


  • Cumulative upstream pressure at BUL-01 is detectable but not severe. BUL-01 integrates the combined influence of agricultural and cattle rangeland use through the valley, the Knockholt Landfill, and the Houston WWTP outfall 500 m upstream. The community retains high EPT richness and HBI (mean 4.14, range 3.91–4.31) is within the “very good — possible slight organic pollution” category, indicating a system experiencing nutrient subsidy rather than severe impairment.

  • Historical context — marine-derived nutrients. The interpretation of “reference” and “enrichment” in this watershed is complicated by the loss of a natural nutrient subsidy. Pacific salmon returning to spawn and die historically delivered substantial marine-derived nitrogen and phosphorus to freshwater systems via carcass decomposition, fueling in-stream and riparian productivity across Pacific Northwest watersheds (Cederholm et al. 1999). The Upper Bulkley once supported significant chinook and sockeye returns (A. S. Gottesfeld and Rabnett 2008), which have declined. Present-day “reference” at BUL-05 therefore reflects a lower nutrient baseline than existed when salmon carcasses were a major annual input, and anthropogenic enrichment at BUL-01 is not simply additive to a pristine baseline — it operates in a system where one nutrient input pathway has diminished while another has grown. The two are biogeochemically distinct: marine-derived nutrients are quickly released from decomposing salmon tissue, pulsed with the timing of salmon runs, and coupled to organic matter and aquatic–terrestrial linkages via scavengers, whereas agricultural, municipal, and industrial inputs are chronic and accompanied by co-contaminants (pharmaceuticals, pesticides, pathogens). The ecological implications of this shift are not resolved by bioindicator data alone.

  • Upstream reaches are in excellent condition, with a caveat. BUL-05’s community structure — high richness, EPT dominance, shredder prevalence, low HBI — represents a high-quality reference condition. However, BUL-05 sits immediately downstream of the McQuarrie Creek confluence. Westcott (2022) placed paired temperature loggers above and below the confluence specifically to quantify McQuarrie’s cold-water thermal influence on the mainstem; Mitchell (1997), conversely, listed McQuarrie among tributaries “highly or severely degraded” by agricultural, municipal, and transportation impacts. These characterizations are not necessarily contradictory — a creek can be thermally beneficial (cool input) and simultaneously nutrient-enriched from upstream land use — but together they leave the net influence of McQuarrie on BUL-05’s benthic community unresolved. The cold, clean thermal signal may create locally favourable conditions at BUL-05 that are not representative of the entire upstream mainstem. Proposed site BUL-06 (Table 8.1), upstream of the McQuarrie confluence, is designed as a Phase 1 test: if the benthic community at BUL-06 is similar to BUL-05, the reference interpretation of BUL-05 strengthens and the current network is sufficient; if BUL-06 instead resembles BUL-04 or BUL-01, the dilution hypothesis becomes the leading explanation and a Phase 2 expansion is triggered — adding a mainstem site immediately above the McQuarrie confluence plus a McQuarrie tributary site to directly measure what the tributary is contributing. Neither Phase 1 outcome fully disambiguates McQuarrie dilution from natural in-stream recovery between BUL-06 and BUL-05, but a dirty BUL-06 is the outcome that makes further investigation actionable and cost-justified.

  • Within-site variability at BUL-04 warrants investigation. The high replicate variability may reflect genuine habitat heterogeneity or could indicate the site is in a transitional zone where small differences in microhabitat produce large differences in community composition.

8.6 Recommendations

  1. Continue annual or biennial monitoring. The three-site design established here captures the key nutrient gradient and can detect temporal changes in community condition with relatively modest sampling effort (nine kick-net samples per campaign).

  2. Establish a routine water quality monitoring program. The biological patterns documented here point to nutrient enrichment, but the existing water quality record is sparse and inconsistent. A structured sampling program — collecting nutrient concentrations (phosphorus, nitrogen species), temperature, pH, dissolved oxygen, and conductivity at a minimum biweekly frequency with additional sampling during high and low flow events — would provide the chemical context needed to interpret benthic community condition and track changes over time. BC Water Quality Guidelines specify a minimum of five samples within 30 days to evaluate long-term chronic guideline compliance (WLRS 2024), making consistent sampling effort essential.

  3. Expand the monitoring network. The current three-site design covers the mainstem nutrient gradient between Houston and McQuarrie Creek but leaves gaps in our understanding of the broader watershed. The Equity Silver Mine straddles the drainage divide between Foxy Creek (flowing north into Maxan Creek) and Bessemer Creek (flowing south into Buck Creek) (Nijman 1996), and both drainage paths warrant monitoring. Existing and proposed monitoring sites are summarized in Table 8.1 and shown in Figure 8.1.


my_caption <<- "Monitoring sites for benthic assessment of the Neexdzii Kwa watershed. Existing sites were sampled in 2025; proposed sites are designed to capture mine drainage, forestry and agricultural influences, point-source impacts, and chinook-bearing reaches."
my_tab_caption(tip_flag = FALSE)
Table 8.1: Monitoring sites for benthic assessment of the Neexdzii Kwa watershed. Existing sites were sampled in 2025; proposed sites are designed to capture mine drainage, forestry and agricultural influences, point-source impacts, and chinook-bearing reaches.
sites_monitoring <- readr::read_csv("data/raw/sites_monitoring.csv", show_col_types = FALSE) |>
  dplyr::filter(!is.na(site_id))

sites_monitoring |>
  dplyr::mutate(reach = paste0(site_id, " \u2014 ", description)) |>
  dplyr::select(Status = status, Stream = stream, Reach = reach, Rationale = rationale) |>
  knitr::kable(booktabs = TRUE) |>
  kableExtra::kable_styling(full_width = TRUE) |>
  kableExtra::row_spec(which(sites_monitoring$status == "Existing"), background = "#f0f0f0")
Status Stream Reach Rationale
Existing Neexdzii Kwa mainstem BUL-05 — Below McQuarrie Creek confluence Cold-water tributary influence. Upstream of Knockholt Landfill
Existing Neexdzii Kwa mainstem BUL-04 — Immediately downstream of Knockholt Bridge (McKilligan Rd) Mid-reach transitional zone; downstream of Knockholt Landfill.
Existing Neexdzii Kwa mainstem BUL-01 — Just upstream of the North Road overpass downstream of Houston adjacent to rest stop. Urban/nutrient influence; matches CABIN historical site MOR37
Proposed Buck Creek BUC-01 — Lowest reach near Bulkley confluence Captures lower reaches of Buck Creek. Utilized chinook spawning and rearing habitat.
Proposed Neexdzii Kwa mainstem BUL-06 — Below Richfield Creek confluence Downstream of legacy concentrate shed site; also a Phase 1 test of whether excellent condition at BUL-05 reflects upstream mainstem condition or localized McQuarrie Creek cold-water input
Proposed Maxan Creek MAX-01 — Upstream of Bulkley Lake. Forestry and agricultural land use influences
Proposed Buck Creek BUC-02 — Buck Ck @ 12 km (CABIN MOR33) Temporal comparison with 2004 CABIN baseline


knitr::include_graphics("fig/map_proposed-sites.png")
Proposed monitoring sites for expanded benthic assessment of the Neexdzii Kwa watershed. Triangles indicate existing 2025 sampling sites; circles show proposed expansion sites. Monitoring site marker colors correspond to the upstream sub-watershed area from which each station receives inputs.

Figure 8.1: Proposed monitoring sites for expanded benthic assessment of the Neexdzii Kwa watershed. Triangles indicate existing 2025 sampling sites; circles show proposed expansion sites. Monitoring site marker colors correspond to the upstream sub-watershed area from which each station receives inputs.