Cyclical, multi-trophic-level responses to a volatile, introduced forage fish: Learning from four decades of food web observation to inform management

Abstract

Species introductions can have significant effects on recipient ecosystems. Anticipating potential ecosystem change in response to introduced species based on historical information can help managers prepare for future conditions. Rainbow Smelt Osmerus mordax have been introduced widely to improve sport fish growth. As intended, Walleye Sander vitreus growth in Horsetooth Reservoir, Colorado increased after Rainbow Smelt introduction, but poor Walleye recruitment occurred as well. Additionally, opossum shrimp Mysis diluviana became absent from both predator diets and intermittent surveys, the dominant Daphnia species in Horsetooth Reservoir shifted and Daphnia densities declined significantly. These patterns were repeated during two different time periods of increased Rainbow Smelt abundance, suggesting that Rainbow Smelt have a strong influence on multiple components of the ecosystem. The repetition of responses to Rainbow Smelt offered the opportunity to evaluate indicators to anticipate potential ecosystem regime shifts that restructure predator–prey dynamics across trophic levels. Three predictors (i.e., high estimated Rainbow Smelt abundance, high catch rates of large Walleye, and low Daphnia densities) were associated with poor Walleye recruitment. Simple indicators like these could inform timely management decisions to take advantage of the benefits Rainbow Smelt offer, while lessening their undesirable effects. For example, management decisions could be made, such as preparing for Walleye egg collections, rearing and stocking of Walleye, increasing availability or quality of Walleye spawning habitat, allowing more protective or liberalized adult Walleye harvest to promote natural recruitment, and limiting Rainbow Smelt access to their spawning habitat.

INTRODUCTION

Understanding and anticipating ecological responses to species introductions is crucial and requires predictive techniques to forecast the potential impacts of invasive species (Dick et al., 2014). There are many examples of predictive models to determine habitat and system suitability for invaders, and where they are likely to spread and establish (e.g., Fournier et al., 2019). However, predicting the deleterious ecological impacts of invasive species establishment and proliferation is less common, especially in aquatic environments (Ricciardi, 2003). Because of the importance of predicting invasive species impacts and anticipating the need for potential management change, Ricciardi (2003) suggested that using empirical data offers opportunity to predict and anticipate potential ecological change from species invasions and proliferation. Additionally, this suggestion was made within the context of empirical data that may lack precision or comprehensiveness, as crude predictions can be informative, and refined further as more contemporary and comprehensive data are collected (Ricciardi, 2003).

Rainbow Smelt Osmerus mordax are an important forage species for a variety of sport fish (Scott & Crossman, 1973). Because of their desirable characteristics, Rainbow Smelt introduction has been used as a management tool to enhance fisheries and encourage sport fish growth (Evans & Loftus, 1987; Mercado-Silva et al., 2006). Thus, Rainbow Smelt have spread from northeastern North America throughout the Great Lakes watershed and the western USA (Mercado-Silva et al., 2006).

Although Rainbow Smelt introductions have been associated with increased growth of sport fish, they have also been associated with undesired food web changes (Mercado-Silva et al., 2007). In some cases, Rainbow Smelt are capable of controlling populations of invertebrate prey species like opossum shrimp Mysis diluviana (Johnson & Goettl, 1999). Further, Galbraith (1967) found that when Rainbow Smelt and Fathead Minnow Pimephales promelas were established in Sporley Lake, Michigan, Daphnia pulex were extirpated and replaced by smaller D. galeata mendotae and D. retrocurva, similar to observations by Reif and Tappa (1966) in Harvey’s Lake, Pennsylvania. Rainbow Smelt also consume and compete with some coregonine species (Evans & Waring, 1987; Hrabik et al., 1998; Loftus & Hulsman, 1986).

Declines in recruitment and abundance of Walleye Sander vitreus have been associated with increasing Rainbow Smelt abundance (Johnson & Goettl, 1999; Mercado-Silva et al., 2007; Schneider & Leach, 1977) which may be a result of predation on (Lepak et al., 2023) and competition with young Walleye for prey resources (Johnson & Goettl, 1999; Lawson & Carpenter, 2014; Mercado-Silva et al., 2007). Further, Rainbow Smelt populations have the potential to exhibit self-regulated, cyclical patterns and alternating year-class dominance, attributable in part to their cannibalistic behavior (He & LaBar, 1994; Henderson & Nepszy, 1989). Thus, by exhibiting strong interactions across multiple trophic levels from the “middle–out” (DeVries and Stein, 1992; also see Stein et al., 1995 for an example involving Gizzard Shad Dorosoma cepedianum), Rainbow Smelt are a fundamental component of many food webs.

In 1983, Rainbow Smelt were introduced into Horsetooth Reservoir, Larimer County, Colorado to enhance growth of Walleye and Smallmouth Bass Micropterus dolomieu (Goettl & Jones, 1984). Historically, Horsetooth Reservoir was managed as a two-tiered fishery with naturally reproducing Walleye (rare in Colorado reservoirs) as a primary component, along with stocked Rainbow Trout Oncorhynchus mykiss and landlocked Sockeye Salmon O. nerka (often called kokanee). The unique, naturally recruiting Walleye had been sustained largely on a forage base of Yellow Perch Perca flavescens, which collapsed, and managers attempted to replace with Rainbow Smelt. Horsetooth Reservoir is the only reservoir in Colorado that supports Rainbow Smelt, and the population has fluctuated widely since their introduction. Because of the importance of the Walleye fishery and reliance on a volatile Rainbow Smelt population, Horsetooth Reservoir is a relatively data rich system with intermittent research and monitoring efforts focused on multiple ecosystem components. Over 40 years of data are available from this system, along with cyclical periods of high densities of Rainbow Smelt during 1987–1996 and since 2013, which provided an opportunity to evaluate long-term, multi-trophic level responses to Rainbow Smelt abundance across two distinct periods. In this paper, we characterize community-level trophic interactions and population dynamics among Rainbow Smelt and other biota, including Walleye, M. diluviana, and Daphnia spp. We used machine learning to identify the best predictors of poor Walleye recruitment based on information from historic and contemporary data sets spanning four decades. Repeated patterns in several ecosystem components along with empirical observations prompted us to hypothesize that indicators linked to Rainbow Smelt abundance would be the best predictors of poor Walleye recruitment. We hoped to identify indices that can inform and prepare fisheries managers for near-term impacts from Rainbow Smelt population dynamics in anticipation of taking actions, such as stocking Walleye fry or fingerlings, increasing availability or quality of Walleye spawning habitat, altering Walleye harvest regulations, and limiting Rainbow Smelt access to their spawning habitat.

Site and historical fisheries management description

Horsetooth Reservoir is a 755-ha impoundment at an elevation of 1,655 m when full (capacity of ∼200 million m³). The reservoir has three distinct basins, is relatively long (∼11 km) and thin (∼1 km) with a mean depth of 25 m and a maximum depth of 70 m (Figure 1). It was completed in 1949, and inflow mainly comes through the 13-mile-long Hansen Feeder Canal from the Colorado–Big Thompson water diversion project.

Following Rainbow Smelt introduction into Horsetooth Reservoir in 1983, Walleye growth increased through 1988 (Jones et al., 1994). Walleye had been reproducing naturally in Horsetooth Reservoir; however, by 1990 Walleye recruitment had precipitously declined. In response, ∼5 and ∼6 million Walleye fry (∼6–7 mm TL) were stocked in mid-April of 1992 and 1993, respectively. Since these stocking efforts appeared unsuccessful, between ∼50,000 and ∼88,000 Walleye fingerlings (∼61–79 mm TL) were stocked annually during July–September 1994–1997. Walleye recruits were subsequently observed in higher numbers and Rainbow Smelt abundance declined. No additional Walleye stocking occurred for the next 22 years (1998–2020). However, in an attempt to offset poor recruitment observed after Rainbow Smelt abundance increased a second time, 3.6 million Walleye fry (∼6 mm TL) were stocked annually during early to mid-April 2021–2023.

While M. diluviana can enter Horsetooth Reservoir through the Hansen Feeder Canal, they were stocked in the reservoir from 1971 to 1974 to enhance the forage base for salmonids (Nesler, 1986). In the 1980s, M. diluviana represented a major prey item for some fish (Nesler, 1986), but were not observed in the reservoir during sampling or in fish stomachs during the two periods when Rainbow Smelt abundance was relatively high—from 1988–2003 and 2010–2015 (Johnson & Goettl, 1999, Jones, 1985a, Silver et al., 2021).

Despite M. diluviana being introduced to enhance salmonid forage, Horsetooth Reservoir does not sustain salmonid natural reproduction, likely due to limited spawning habitat and whirling disease. Therefore, ∼350,000 kokanee fry were stocked annually as well as ∼40,000 catchable (> 150 mm TL) and 225,000 sub-catchable (< 150 mm TL) Rainbow Trout were stocked prior to 2000 to sustain a sport fishery. After 2000, kokanee stocking was ended, and Rainbow Trout stocking was reduced to an average of ∼13,000 catchable trout and 22,000 sub-catchable trout annually. Additionally, based on previous research in Colorado reservoirs, it is likely that adult Walleye consume most of the stocked Rainbow Trout (Lepak et al., 2012; Stacy & Lepak, 2012). However, little information is available on salmonid distribution and abundance in Horsetooth Reservoir, since routine sampling surveys were designed to capture Walleye and other nonsalmonid sport fish.

Horsetooth Reservoir fluctuates significantly (on the order of 10–30 m annually; see Supplemental Information; Figure S1) since it is used as a municipal and agricultural water source. This limits the littoral habitat available for aquatic organisms. For example, at the end of 2000, the reservoir was drawn down substantially to < 5% (∼8 million m³) capacity for maintenance, and at times during the winter was over 40 m below full pool (Northern Colorado Water Conservancy District, unpublished data). In the early spring of 2001–2003, Horsetooth Reservoir was at < 10% capacity, limiting Rainbow Smelt and Walleye access to known spawning habitat in the reservoir inlet. However, there were indications that Rainbow Smelt had declined in the years leading up to the drawdown (Johnson & Goettl et al., 1999) and their decline continued through 2010. After 2010, the Rainbow Smelt population increased again, and the pattern of low Walleye recruitment, low M. diluviana density, and low Daphnia density, and dominance by D. galeata was ­repeated.

METHODS

Walleye diet and stable isotope analyses

Gill netting in late April through May and boat electrofishing in July and August surveys are done intermittently in Horsetooth Reservoir (see more detail in section “Walleye recruitment and abundance indices” below). Individual Walleye were retained from these sampling events to examine diet and growth.

Walleye diets from 1983–1992 and 1994–1996 were based on stomach contents of collected fish. Stomach contents were removed and analyzed under a dissecting microscope to identify each prey item to taxon, and mean percent volume was calculated for each prey item to obtain a snapshot of what fish had recently eaten (Jones, 1985a). In 2008, carbon and nitrogen stable isotope values of Walleye and their prey were used to obtain a time-integrated estimate of diet (Johnson et al., 2015; Post, 2002). These data were analyzed using Bayesian mixing models. In 2013 and 2017–2019, these data were analyzed with a model written in R that is designed to solve mixing equations for stable isotopic data within a Bayesian framework (stable isotope mixing models in R [SIMMR] 0.4.5; Parnell, 2021). Trophic fractionation (i.e., differences in carbon and nitrogen stable isotopes between prey and predator) and error for δ¹³C (0.4 and 1.3, respectively) and δ¹⁵N (3.4 and 1.0) were set to values established by Post (2002). Values for burn-in (1,000 separate simulations) and the number of iterations (each simulation run 10,000 times) were the SIMMR default settings, and additional iterations were unnecessary, since Gelma–Rubin convergence diagnostics never exceeded 1.01 for any model.

In 2013, 30 Walleye were collected for analyses. Prey species targeted for collection had to be abundant, of edible size, and common in Walleye diets in the past. These prey items were collected opportunistically in Horsetooth Reservoir using beach seines, various nets, minnow trapping, and boat electrofishing. Prey items included in analyses were Rainbow Smelt (n = 4), Decapoda (n = 6), hatchery-reared Rainbow Trout (n = 9), dipterans (n = 5 composite monthly samples of multiple individuals collected from May to October), large zooplankton (n = 3 composite samples of multiple individuals representing bulk zooplankton collected with a 500 micron mesh net in June, July, and September), Gizzard Shad (n = 3), and Yellow Perch (n = 3). Prey items were placed into the following categories: Rainbow Smelt, Decapoda (crayfish), salmonidae (Rainbow Trout), other invertebrates (dipterans and zooplankton), and other fish (Gizzard Shad and Yellow Perch).

In 2017, 2018, and 2019, 42, 14, and 42 Walleye were collected for analysis, respectively. During these years, prey species targeted for collection had to be abundant, of edible size, and common in Walleye diets in the past. Cursory diet analyses indicated that Walleye were mainly focused on consuming Rainbow Smelt. Thus, Rainbow Smelt were targeted with fine mesh gill nets all three years, and other prey species were collected opportunistically with netting, trapping, and boat electrofishing efforts. Prey items collected for analyses were Rainbow Smelt (n = 213 across all three years), Decapoda (n = 6), large zooplankton (n = 6 composite samples of multiple individuals representing bulk zooplankton collected monthly with a 500-micron mesh net from June through November 2019 with values averaged across samples from each of the three reservoir basins), and Gizzard Shad (n = 36 from 2017 and 2018). These prey items were placed into the following categories: Rainbow Smelt, Decapoda, other invertebrates (zooplankton), and other fish (Gizzard Shad). Because samples from 2017–2019 were collected during a relatively short time window, and values for Rainbow Smelt stayed relatively consistent during this time, prey items from all three years were used for mixture modeling of Walleye.

Stable isotope values were measured using a Thermo Delta V isotope ratio mass spectrometer interfaced to a NC2500 elemental analyzer at the Cornell Isotope Laboratory (Ithaca, New York). In-house standards were analyzed every 10 samples to ensure precision and accuracy. The overall standard deviation of these samples fell below 0.20‰ for δ¹³C and δ¹⁵N. A methionine standard was used to quantify the ability of the instrument to measure across a gradient of amplitude intensities. Values for δ¹³C and δ¹⁵N were corrected using primary references of Vienna Pee Dee Belemnite (a Cretaceous marine fossil, Belemnitella americana, from the Pee Dee Formation in South Carolina used as a standard for carbon-13), and atmospheric air. To avoid potential bias from differing lipid concentrations among samples and species, corrections for lipid ­content from Post et al., (2007) were applied to δ¹³C values:

where C:N is the carbon-to-nitrogen ratio.

Rainbow Smelt diet and stable isotope analyses

Rainbow Smelt were collected using fine and large mesh gill nets, trawls, boat electrofishing, and beach seines throughout the study period. Rainbow Smelt diets were analyzed using stomach contents of fish collected in 1984 (Jones, 1985a), 1987 (Thomas, 1989), 1989 to 1992 (Goettl, 1990, 1991, 1992, 1993), 1994 (Goettl & Johnson, 1995) and 1995 (Goettl & Johnson, 1996), or stable carbon and nitrogen isotope values of predator and prey coupled with a Bayesian mixing model (SIMMR 0.4.5) for samples collected during 2017–2019. For both volumetric diet composition and stable carbon and nitrogen isotope analyses, Rainbow Smelt prey were grouped into four categories: M. diluviana, Rainbow Smelt, zooplankton (e.g., daphnids, copepods), and other invertebrate prey (e.g., dipterans, amphipods).

For stomach content analyses, Rainbow Smelt ≥ 100 mm TL were commonly analyzed separately due to their potential to consume M. diluviana, Walleye larvae, and smaller Rainbow Smelt. In 1984, 1987, 1989, and 1991, only data from Rainbow Smelt ≥ 100 mm TL were used. In 1984, two size-classes of 100–150 mm TL (n = 3) and > 150 mm TL (n = 2) were considered. In 1987, three size-classes of 100–125 mm TL (n = 16), > 125–150 mm TL (n = 21), and > 150 mm TL (n = 51) were used. Diet composition was weighted by sample sizes within each respective size-class. In 1989, all Rainbow Smelt > 100 mm TL (n = 90) were used together, and the same was done in 1991 (n = 99). In 1990, 1992, 1994, and 1995, data were collected from Rainbow Smelt sampled to all sizes (n = 94, 212, 435, and 131, respectively). Stomach contents were removed and examined under a dissecting microscope, identified to taxon, and quantified as mean percent volume by prey item (Jones, 1985a).

For stable isotope analyses in 2017–2019, we classified Rainbow Smelt as either predator (≥ 130 mm TL) or prey (< 130 mm TL) since we noted a change in isotope values, which indicated a diet shift, when Rainbow Smelt reached 130 mm TL. There were 20, 9, and 73 Rainbow Smelt ≥ 130 mm TL in 2017, 2018, and 2019 samples, respectively. Other items measured and included in the stable isotope mixing models were zooplankton (as described in the Walleye diet analysis), and dipterans, which were regularly found in Rainbow Smelt stomachs, but not collected during these years. Thus, we included data from dipterans collected in 2013 (described in the Walleye diet analysis) despite disparate collection time periods. In 2019, we collected Rainbow Smelt (n = 500) using suspended fine mesh gill nets from May to November. All stomach contents from these fish were inspected for the presence of M. diluviana and larval and juvenile Walleye. In all cases (volumetric diet composition and composition established with stable carbon and nitrogen isotope analyses) Rainbow Smelt prey categories were collapsed into M. diluviana, Rainbow Smelt, zooplankton (e.g., daphnids, copepods), and other invertebrate prey (e.g., dipterans, amphipods).

Walleye recruitment and abundance indices

Gill netting (late April–May) surveys have occurred most years beginning in 1983 during May in Horsetooth Reservoir (Jones et al., 1994). Generally, a set of 20 standardized monofilament gill nets are used to assess the adult Walleye population. Nets are ∼46 m long and ∼1.8 m in depth, with 6 graduated panels ranging from 13–76-mm bar-mesh increasing in 6–13-mm increments and set overnight in the same location perpendicular to shore. Net soak times were often estimated to the nearest hour, reported to the nearest minute, or reported as “overnight,” which limited us to using this metric (“overnight”) when combining these data. Thus, Walleye catch per unit effort (CPUE) during spring gill netting surveys was expressed as Walleye catch per overnight net set. Catches from each of the 20 nets combined were used to calculate CPUE for that year. To remain consistent with previous research efforts (e.g., Johnson & Goettl, 1999), two different size-classes of Walleye were used to characterize the population by representing reproductive adult Walleye and juvenile Walleye recruiting to the sampling gear: 150–300 mm TL (predominately age-1 Walleye recruits based on our unpublished data), and ≥ 451 mm TL (predominantly adult reproductive Walleye). While it would have been preferable to also assess young-of-the-year Walleye, the gill net mesh size used was too large to do this. Therefore, we opted to use the most consistent long-term source of Walleye catch data to develop indices and predictors.

Walleye growth

Historically, Walleye growth comparisons were standardized to mean TL at age 3 (Goettl & Jones, 1984, 1986), so we used this metric for the entire time-series. Walleye TL at age 3 in 1991 and earlier were back-calculated from scales collected in 1983 (n = 204), 1984 (n = 210), 1985 (n = 94), 1986 (n = 93), 1988 (n = 68), 1989 (n = 97), 1990 (n = 138), 1991 (n = 94), and 1992 (n = 187; Goettl, 1990, 1991, 1992, 1993; Goettl & Jones, 1984, 1986; Goettl & Thomas, 1987; Jones, 1985b). Walleye TL at age 3 in 2001 and later were back-calculated from sectioned otoliths (Farrell et al., 2021) collected in 2012 (n = 40), 2013 (n = 119), 2019 (n = 131), and 2021 (n = 75). Mean Walleye length at age 3 from scales and otoliths were computed each year when ≥ 2 individuals were represented. Standard error was calculated from 2001–2016, but only mean ages were available for earlier data, and these were weighted by sample size to estimate yearling means.

Rainbow Smelt abundance indices

When available, hydroacoustic surveys were used to estimate Rainbow Smelt abundance. During 1993–1996, surveys were conducted along longitudinal transects down the center of the north–south axis using a BioSonics 420-kHz dual-beam echosounder and analyzed using echo integration with in situ trawl catch data to determine mean Rainbow Smelt body size and target strength (-50 dB; Johnson & Goettl, 1999). Due to lower Rainbow Smelt abundance during 1998–2013, surveys during those years were conducted along longitudinal transects down the center of the north–south axis, but with a 200-kHz Hydroacoustic Technology model 243 split beam digital echosounder and real-time target tracking algorithm. Individual fish targets were identified based on a minimum of four echoes, and target strengths corresponding to individuals > 5 cm TL and < 20 cm TL were based on fine mesh gill net catches (our unpublished data) and a target strength–length relationship for Rainbow Smelt (Rudstam et al., 2003). Rainbow Smelt were largely observed near the thermocline and in the epilimnion during sampling based on hydroacoustics and vertical gill netting. Targets below 30 m, which were usually fish > 20 cm TL, were not used in analyses. During 2017–2023, surveys were conducted along zig-zag transects down the center of the north–south axis of the reservoir with a 200-kHz Hydroacoustic Technology model 243 split beam digital echosounder (the same equipment as 1998–2013), and echo integration was used instead of target tracking to enumerate Rainbow Smelt, and in situ target strengths were used to determine Rainbow Smelt density (Johnson & Goettl, 1999). We note that gear, personnel, and methodology were ­inconsistent between historic and contemporary hydroacoustic surveys, but other data are unavailable from this system during the time periods of interest.

Prior to initiation of hydroacoustics surveys in 1993, Rainbow Smelt abundance was documented using midwater trawling surveys. Trawl surveys were conducted at night in a stepped-oblique approach from the surface to a depth of 18 m using a net with a 6 × 6-m opening and graduated mesh from 203 mm (stretch) at the open end to 13 mm at the cod end (Kirn & LaBar, 1991). The relationship between estimated Rainbow Smelt densities from these surveys and hydroacoustic surveys in years these surveys overlapped are described in Johnson & Goettl (1999). We also used routine gill netting (late April–May) and boat electrofishing in July–August (Jones et al., 1994) to develop an annual index of Rainbow Smelt based on data from different gears. The coarse index was represented by the total Rainbow Smelt catch across these surveys divided by the number of surveys that occurred in a given year. This index was calculated because we were interested whether Rainbow Smelt catches during surveys not designed to target Rainbow Smelt would be still be able to identify times of high Rainbow Smelt abundance. This approach would allow these surveys to replace hydroacoustic and fine mesh gill netting sampling, which would save time, effort, and money.

Mysis diluviana density

Estimates of M. diluviana density (number/m²) were obtained at 10 locations using a 1-m diameter, 500-μm mesh net that was towed vertically through the water column at night in 1999, 2003–2006, 2008, 2012, 2013, 2015–2018, and 2020–2023 in strata representative (proportionally by depth) of regions in Horsetooth Reservoir that support M. diluviana. Specifically, strata included 10–20 m (n = 1), > 20–30 m (n = 1), 30–40 m (n = 4), 40–50 m (n = 3), and 50–60 m (n = 1) depths depending on water conditions (Silver et al., 2021). Benthic trawl data from 1981 were used to develop an early, conservative M. diluviana density estimate (Nesler, 1986).

Daphnia density and dominant species

Horsetooth Reservoir is similar to some other Colorado Reservoirs in that the highest densities of zooplankton generally occur in surface waters (< 10-m depth; e.g., Johson and Goettl, 1999; Nelson, 1971). During 1987–1993, macrozooplankton sampling in Horsetooth Reservoir was conducted April–September in mid-morning with a metered Clarke–Bumpus net fished in the top 5 m of the water column (Johnson & Goettl, 1999). During 1994–1996, 2002–2023, sampling was also conducted April–September in mid-mornings, but with a 153-μm Wisconsin net fished in the top 5 m or 10 m of the water column. Samplings during these years were done by either Colorado State University, the Bureau of Reclamation and Northern Colorado Water Conservancy District, or Colorado Parks and Wildlife personnel. Since some methodology changed over time, we used the sampling day from each year with the highest estimated mean density of Daphnia (individuals/L) as an index of macrozooplankton abundance. If sampling was done in more than one of the reservoir’s three basins during the same year, we averaged data from all basins sampled for each sampling day. We also calculated the percent of large-bodied D. pulex/pulicaria within each year’s index sample.

Water storage index

Since low water levels can negatively influence access to both Rainbow Smelt and Walleye spawning areas in Horsetooth Reservoir, water storage during spawning periods was considered as a factor that could influence relationships between Walleye recruitment and Rainbow Smelt. During the 2001–2003 drawdowns in Horsetooth Reservoir, limited access to known spawning habitat for Rainbow Smelt and Walleye in the reservoir inlet was observed (K. Kehmeier, Colorado Parks and Wildlife, personal communication). Sampling, observation, and egg-take operations suggest Rainbow Smelt and Walleye begin to spawn in late March. So, using data provided by the Northern Colorado Water Conservancy District, we calculated annual mean daily March water storage as an indicator of water ­conditions.

Machine learning prediction

We compiled disparate abiotic and biotic data from Horsetooth Reservoir and developed annual indices to evaluate their importance as predictors of Walleye recruitment failure/success. We applied a machine learning approach (Random Forest; RF) to evaluate the importance and performance of seven predictors of Walleye recruitment:

1. Gill net catch rates of large (≥ 451 mm TL) Walleye (LWAL; Walleye overnight/net).

2. Estimated number of Rainbow Smelt from hydroacoustic surveys (SONAR).

3. Rainbow Smelt encountered during routine sampling efforts (RSM; number/survey).

4. Estimated M. diluviana density (MYSIS; M. diluviana/m²).

5. Estimated peak Daphnia density (DAPHd; Daphnia/L).

6. The percent of peak Daphnia sample composed of larger-bodied D. pulex/pulicaria (DAPHc; %).

7. Reservoir water storage (H2O; km³).

A one-year offset was applied to the CPUE of Walleye 150–300 mm TL during routine spring gill netting to coincide with data collected the previous year when that cohort of Walleye hatched (to reflect hatching conditions). Missing data were predicted using linear interpolation (detailed methodology available in Supplemental Information; “Linear ­interpolation”).

A machine learning approach was selected because many of the predictors were highly correlated and could have complex interactions that precluded the use of traditional ­statistical approaches. Significant (p < 0.01) Pearson correlation coefficients (r) were found between LWAL and SONAR (r = 0.72); LWAL and MYSIS (r = -0.41); LWAL and DAPHd (r = -0.44); and LWAL and DAPHc (r = -0.57); LWAL and H2O (r = 0.44); SONAR and MYSIS (r = -0.48); SONAR and DAPHd (r = -0.47); SONAR and DAPHc (r = -0.58); SONAR and H2O (r = 0.43). To perform the RF analysis, we used the ‘randomForest’ package in R (R Development Core Team 2011; version 4.2.2). Accuracy and error rate calculations for each observation using out-of-bag predictions were based on 2,000 regression trees (i.e., predicting data withheld from each tree). Data predictions that were not used to evaluate fit were considered a form of cross-validation. Variable importance was assessed by comparing the increase in mean squared prediction error and node purity associated with each individual covariate (Cutler et al., 2007; Gini, 1912; Liaw & Wiener, 2002). Increases in mean square error are reflected as increases in Walleye recruitment estimate error, while increases in node purity were used to assess variable importance. At each node split in a tree, there is a resulting decrease in the Gini index (a measure of the discrepancy between observations and predictions; Gini, 1912). The sum of these decreases throughout the RF for a given variable, normalized by the number of trees was used as an indicator of variable importance (node purity). Higher node purity indicated that a variable was more important for Walleye recruitment prediction when compared to another.

Based on results, we opted to develop partial dependence plots of Walleye 150–300 mm TL CPUE as a function of the observed values for the three best predictors (SONAR, LWAL, and DAPHd). This approach provides predictions associated with a single variable at observed points while averaging, and holding constant, the effects of other predictors to better visualize the effect of the variable of interest (pdp Package in R; Greenwell, 2017).

RESULTS

Trophic and predator–prey interactions

In 1984, the stomach of one 255-mm TL Walleye collected during the routine gill net survey contained Rainbow Smelt (Jones, 1985a). The presence of Rainbow Smelt in Walleye diets from routine surveys increased thereafter and were a significant component of Walleye diets from 1988–1996, along with invertebrates and salmonids (e.g., Rainbow Trout), which were stocked more frequently prior to 2000 (Johnson & Goettl, 1999; Jones et al., 1994). Walleye largely consumed crayfish when Rainbow Smelt were less abundant (Johnson et al., 2015; K. Kehmeier, personal communication) and salmonid stocking was reduced. This information, coupled with stable carbon and nitrogen isotope data in 2008, 2013, and 2017–2019, indicated that Walleye consumed a mix of prey sources in the reservoir in 2008, but transitioned to primarily a Rainbow Smelt diet by 2017 (Figures 2, 3). Other invertebrates were generally 5% or less of Walleye diets during 2017–2019 based on the mixing model approach (Figures 2, 3).

When Rainbow Smelt were first introduced into Horsetooth Reservoir, M. diluviana composed most of their diet. Mysis diluviana disappeared from Rainbow Smelt diets by 1988 (Johnson & Goettl, 1999), and were replaced by Daphnia and dipterans, as well as an increasing occurrence of cannibalism in the 1990s after first being observed in 1987 (Figure 4). Only zooplankton were observed in Rainbow Smelt diets in 1994 and 1995 (Goettl & Johnson, 1995, 1996), but in 2017–2019, stable carbon and nitrogen values indicated that Rainbow Smelt ≥ 130 mm TL were getting large portions of their energy from cannibalism (Figures 3, 4). During more contemporary stable isotope analyses, M. diluviana were not readily available for collection, so they could not be included in the SIMMR analyses. Further, no M. diluviana or larval Walleye were observed in 500 Rainbow Smelt stomachs from fish collected across seasons in 2019 from May to November, and most prey biomass was Rainbow Smelt. Analyses of mitochondrial DNA from stomach contents of Rainbow Smelt collected near the inlet of Horsetooth Reservoir in 2022 indicated they were consuming some Walleye eggs or larva (Lepak et al., 2023). Rainbow Smelt diets linked them to organisms within multiple trophic levels from the middle–out, corroborating previous findings and allowing for development of a predictive framework (DeVries and Stein, 1992; Johnson & Goettl, 1999; Stein et al., 1995).

Walleye recruitment and abundance indices

Based on standardized gill net catches, Walleye 150–300 mm TL were rarely encountered during periods of relatively high Rainbow Smelt abundance. Catch rates of these small Walleye were generally below 0.5 fish/net during these periods, indicating poor recruitment (Figure 5A). Concomitant with Walleye fingerling stocking during 1994–1997 to supplement natural reproduction, Rainbow Smelt abundance declined and catch rates of Walleye 150–300 mm TL increased to > 0.5 fish/net (Figure 5A). Catch rates occasionally approached or exceeded 2.0 fish/net until 2012, when the Rainbow Smelt population rebounded and catch rates of small Walleye declined again (Figure 5A).

Catch rates of large Walleye (≥ 451 mm) in standardized gill net samples showed increases during periods of elevated Rainbow Smelt abundance (Figure 5B). The pattern of increasing catches of large Walleye was observed in the 1990s when Rainbow Smelt were abundant (Johnson & Goettl, 1999), and then again beginning around 2010. This repetition was expected, as consumption of Rainbow Smelt forage previously produced relatively large Walleye in Horsetooth Reservoir (Johnson & Goettl, 1999; Jones et al., 1994), and likely increased the number of large fish available for capture.

Walleye growth

Walleye growth responded positively to increases in Rainbow Smelt abundance in the late 1980s to mid-1990s (Figure 5C). Walleye growth after the introduction of Rainbow Smelt appeared to remain high (an increase in length by 50% at age 3) compared to observations prior to Rainbow Smelt introduction (Jones et al., 1994). Notably, poor Walleye recruitment created large cohort gaps (e.g., periods after 1989 and 2013) and few Walleye were captured from cohorts following the first Rainbow Smelt population increase and leading up to the second (Figure 5). These gaps and low sample sizes are informative because they indicate large declines in Walleye recruitment and year-class failure lagging approximately 5 years after Rainbow Smelt abundance increased.

Rainbow Smelt abundance indices

Hydroacoustic and trawl surveys indicated there were millions of Rainbow Smelt soon after their introduction in Horsetooth Reservoir through the mid-1990s (Figure 5D). Rainbow Smelt estimates peaked at almost 10 million individuals in 1994 during their initial expansion period but declined by 2000 to some of the lowest abundances estimated with hydroacoustics, and these observations were 2 and 3 years before drawdowns. After Horsetooth Reservoir was drawn down during 2001–2003, Rainbow Smelt were not observed in samples until 2010, when a single Rainbow Smelt was captured during a gill net survey. Based on hydroacoustic data, Rainbow Smelt abundance began increasing again after 2010, and by 2018, they achieved numbers similar to those observed during their first expansion. Rainbow Smelt catches during routine spring and summer surveys also indicated there were two pulses of increased abundance (Figure 5E).

Mysis diluviana density

Mysis diluviana presence in Horsetooth Reservoir was first confirmed in 1981 from benthic trawl catches (Nesler, 1986). The presence of M. diluviana was also apparent in fish diets (e.g., Rainbow Smelt ≥ 100 mm TL) until 1988 (Jones, 1985a; Thomas, 1989). Subsequently, very few M. diluviana were captured until after the 2004 drawdown and refilling of Horsetooth Reservoir. For the next 10 years, during a period of relatively low Rainbow Smelt abundance, M. diluviana appeared more frequently during sampling, and surveys indicated their abundance had increased relative to observations around 2000 (Figure 5F). However, M. diluviana dropped below detection limits again by 2020 following the resurgence of Rainbow Smelt abundance.

Daphnia density and dominant species

Peak density of Daphnia did not exceed 10 individuals/L ­during 1989–1994. During this period, small-bodied ­D. ­galeata mendotae was the dominant species. This pattern of low Daphnia/L dominated by D. galeata mendotae occurred again in Horsetooth Reservoir during 2012–2023 (Figure 5G). Both these time periods corresponded to relatively high Rainbow Smelt abundance. Daphnia pulex/pulicaria are generally considered more desirable prey items because they tend to be larger than D. galeata mendotae. Thus, it appears that when Rainbow Smelt abundance is high, Daphnia abundance and species composition changes as described here and previously during 1989–1994 (Johnson & Goettl, 1999). In contrast, when Rainbow Smelt abundance is relatively low, higher densities of Daphnia (≥ 10 individuals/L) usually dominated by D. pulex/pulicaria are observed (Figure 5G).

Water storage index

Horsetooth Reservoir water storage fluctuated significantly over the course of the study (Supplemental Information; Figure S1). During 2001–2003, the reservoir was drawn down significantly during Rainbow Smelt and Walleye spawning periods (Figure 5H). Rainbow Smelt had already shown signs of decline prior to this drawdown in 1998 and 1999 (Figure 5D). Walleye recruitment indices following the drawdown remained relatively high for several years, while Rainbow Smelt abundance indices remained low (Figure 5A, 5D).

Machine learning prediction

Rainbow Smelt abundance estimates from hydroacoustic surveys (SONAR), the number of large Walleye captured during routine spring gill net surveys (LWAL), and peak density of Daphnia (DAPHd), were the three most important predictors of the Walleye recruitment index. These three predictors had higher importance metrics compared to the other four predictors examined (Figure 6). Overall, the RF approach performed relatively well with a mean square of residuals equal to 0.39 and explaining 55% of variance in the observed data (Figure 5A). Typical of RF, the model tended to over predict low Walleye recruitment index values and under predict high values (Figure 7). Based on partial dependence plots, we found that poor Walleye recruitment was correlated with Rainbow Smelt abundance estimates greater than 3 million (Figure 8A), large Walleye catch rates greater than 1 fish/net (Figure 8B), and Daphnia densities less than 10 individuals/L (Figure 8C). Conversely, high Walleye recruitment was correlated with Rainbow Smelt abundance estimates less than 3 million, large Walleye catch rates greater than 1 fish/net, and Daphnia densities of at least 10 individuals/L.

DISCUSSION

Predicting deleterious impacts from invasive species to anticipate potential ecological change is important from a management perspective. For example, extensive planning and resources may be necessary to respond to the introduction and establishment of invasive species (e.g., Dreissenid mussels in the Great Lakes region; Nalepa & Schloesser, 2013). Thus, forecasting ecological change before it occurs allows for time and resource allocation to prepare for and respond to undesirable ecological and economic impacts from invasive species. Predicting poor recruitment failure is of particular importance to Horsetooth Reservoir managers because Walleye eggs are not regularly collected and propagated for stocking when Walleye are reproducing naturally (K. Kehmeier, personal communication). If indicators the previous year suggest that recruitment of Walleye the following spring will be limited, actions can be taken to collect and raise Walleye eggs to the fry or fingerling stage for stocking. Allocating resources and preparing for egg collection and propagation is one management response to poor Walleye recruitment. However, managers may also consider other actions like increasing or improving Walleye spawning habitat (which could require reservoir operation changes), altering harvest regulations to protect spawning adult Walleye, encouraging harvest of large Walleye when they cannibalize sublegal juveniles, or limiting Rainbow Smelt access to their spawning habitat (i.e., blocking the reservoir inlet during spawning) to reduce their abundance. All of these actions require resources and forethought to plan and implement. Thus, anticipating ecological change can be important for balancing valuable fisheries and management actions in the context of invasive species and their deleterious impacts.

Rainbow Smelt exhibited strong middle–out influence on the Horsetooth Reservoir food web since their introduction in 1983. Significant responses by organisms at trophic levels above and below Rainbow Smelt were noted when the Rainbow Smelt population increased, and reciprocal observations were made when the population declined (Johnson & Goettl, 1999; Jones et al., 1994). We observed the influence of increasing Rainbow Smelt abundance on Horsetooth Reservoir biota during two separate time periods, providing an opportunity to evaluate biotic responses iteratively to corroborate initial ­observations. Our observations demonstrated that Walleye, M. diluviana, and Daphnia exhibited repeated responses to high Rainbow Smelt abundance, suggesting Rainbow Smelt play a dominant role in driving food web interactions and population dynamics across trophic levels in Horsetooth Reservoir.

Rainbow Smelt have been associated with Walleye recruitment failure in other systems (Evans & Loftus, 1987; Mercado-Silva et al., 2007). This undesirable outcome also occurred in Horsetooth Reservoir, and Walleye recruitment essentially ceased when Rainbow Smelt abundance was high (Johnson & Goettl, 1999). Thus, high Rainbow Smelt abundance appears to improve Walleye growth, but comes at the expense of Walleye recruitment. Rainbow Smelt consumption of larval Walleye near the inlet of Horsetooth Reservoir was observed (Lepak et al., 2023) and is likely reducing Walleye recruitment directly. However, this finding does not discount other potential causes of poor recruitment. For example, the presence of relatively high densities of large Walleye coincided with periods of poor Walleye recruitment, and this could be attributed to cannibalism or competition among high densities of larvae, creating a Ricker-type stock–recruitment relationship (Ricker, 1954). However, despite observations of Walleye cannibalism in other systems (e.g., Forney, 1980; Hansen et al., 1998) information on piscivory by Horsetooth Reservoir Walleye indicates they focused largely on Rainbow Smelt when they were abundant based on diet information. Conversely, when Walleye larvae are at their lowest density (during periods of poor recruitment), Rainbow Smelt are likely available in high densities as alternative Walleye forage. Therefore, it appears that Walleye cannibalism is a less important driver of poor Walleye recruitment when compared with competition or consumption of Walleye by Rainbow Smelt, and the presence of large Walleye is an indicator of poor recruitment rather than the cause. Little is known about larval Walleye density or competition for prey in Horsetooth Reservoir, though larval Walleye densities were low (∼0.5–3.5 larvae/1,000m³) during ichthyoplankton sampling near the inlet in 2022 (Lepak et al., 2023). Interestingly, RF model results predicted exceptionally poor Walleye recruitment when peak Daphnia densities were below 10 individuals/L, and zooplankton was dominated by D. galeata mendotae rather than D. pulex and pulicaria. Walleye recruitment indices were below 0.25 fish/net in these cases, which was similar to observations during periods of what was considered Walleye recruitment “failure” historically.

Intraspecific interactions between Rainbow Smelt include competition for prey resources and cannibalism of small Rainbow Smelt by larger individuals. Rainbow Smelt cannibalism has been observed widely (e.g., Evans & Loftus, 1987; Jones et al., 1994), and may contribute to self-regulation of population density (He & LaBar, 1994; Henderson & Nepszy 1989). It was evident based on diet data, especially at smaller sizes, that Rainbow Smelt were effective planktivores in Horsetooth Reservoir and likely compete strongly with other planktivorous species and life stages of fish. A distinct enrichment was noted in Rainbow Smelt stable nitrogen isotope values at 130 mm TL, indicating a shift to cannibalism within Horsetooth Reservoir. Thus, large Rainbow Smelt ­affect recruitment within their own population by functioning as predators and competitors of smaller Rainbow Smelt.

Mysis diluviana also responded strongly to Rainbow Smelt introduction in Horsetooth Reservoir. Rainbow Smelt consume Daphnia and can compete for zooplankton resources with M. diluviana, but based on diet information, Rainbow Smelt also consume M. diluviana directly. However, M. diluviana can be highly herbivorous or take advantage of detrital resources when zooplankton prey are scarce (Hansen et al., 2023). Thus, M. diluviana are likely affected more through Rainbow Smelt predation than resource competition. These factors, in addition to the intermittent presence of a low dissolved oxygen layer near the thermocline (Silver et al., 2021) appear to drive the M. diluviana population within Horsetooth Reservoir below detection levels when Rainbow Smelt abundance is high. We are unaware of other reports of such a strong response by M. diluviana in connection with a fish species (e.g., Bruel et al., 2021; Pothoven et al., 2009).

Macrozooplankton populations in Horsetooth Reservoir were restructured by Rainbow Smelt. In particular, the peak density of Daphnia in surface waters dropped below 10 individuals/L, and was accompanied by a complete shift in Daphnia species from D. pulex and pulicaria to generally smaller-bodied D. galeata mendotae. Similar changes in zooplankton communities and their densities in response to Rainbow Smelt have been reported by other investigators (e.g., Galbraith, 1967; Reif & Tappa, 1966). It also appears Daphnia density peaks are restricted to periods of peak thermal stratification in mid-to-late summer, which offers Daphnia some thermal refuge from Rainbow Smelt and M. diluviana (Johnson & Goettl, 1999). However, these periods occur after the time larval Walleye need to transition to consumption of larger, preferred Daphnia prey (Johnson & Goettl, 1999). Thus, in Horsetooth Reservoir, macrozooplankton restructuring may be creating conditions where organisms like larval Walleye do not have access to prey resources (May et al., 2021).

Water storage was a poor predictor of Walleye recruitment. While it is possible that low water levels could have disadvantaged Rainbow Smelt disproportionately to Walleye, this was not examined during the reservoir drawdown. Rainbow Smelt did decline prior to the reservoir drawdown, and it is possible declines were perpetuated by low water. It is also possible that the longer generation time in Walleye (Scott & Crossman, 1973), or a lower reliance of Rainbow Smelt on inlet spawning habitat could have allowed Walleye to persist during low water periods and subsequently reproduce and recruit successfully. Walleye recruitment indices in 2001 and 2004 were the highest observed in over a decade during and following the reservoir drawdown. Walleye and Rainbow Smelt recruitment in relation to future similar drawdowns should be studied to help to determine whether similar effects occur.

Unintended consequences occurred as a result of Rainbow Smelt stocking in Horsetooth Reservoir. In particular, periods of elevated Rainbow Smelt abundance led to poor Walleye recruitment, and the need to adapt management strategies. Using food web indicators assessed during routine monitoring, managers can preemptively prepare and allocate resources in response to fluctuations in Rainbow Smelt ­populations.

The predictive framework we developed using an RF model relied on a variety of indices of different quality and precision. The model identified the best three predictors of Walleye recruitment as the estimated number of Rainbow Smelt in Horsetooth Reservoir during hydroacoustic surveys (SONAR), the CPUE of large Walleye (LWAL), and the peak density of Daphnia (DAPHd), and some of these predictors (specifically SONAR and DAPHd) changed prior to observing extended periods of poor Walleye recruitment. We acknowledge conditions in Horsetooth Reservoir are unique, but managing in the presence of introduced species is common around the world (Ricciardi & Simberloff, 2009). Though predictions can be refined with more or higher quality data, this type of approach could provide managers of other species and systems with the ability to anticipate and prepare for ecological change, particularly if responses to introduced species seem to follow predictable patterns or cycles.

SUPPLEMENTARY MATERIAL

Supplementary material is available at Fisheries online.

Figure S1. Horsetooth Reservoir storage. The black line indicates daily water storage (km3) in Horsetooth Reservoir through time. Drawdown conditions are notable from 2001-2003. The plot begins 1 January 1980.

DATA AVAILABILITY

Data are available upon reasonable request.

ETHICS STATEMENT

Sampling was approved by Colorado Parks and Wildlife, and care and use of experimental animals was done in compliance with the guidelines and policies approved by the Colorado Parks and Wildlife scientific collection permit DOW087.

FUNDING

Colorado Parks and Wildlife and Colorado State University personnel provided the primary resources and support for this research.

ACKNOWLEDGMENTS

We thank Miranda Moll and Travis Hackett for their support during laboratory analyses. We thank Grant Wilcox for map production. We also thank Colorado Parks and Wildlife, and Colorado State University personnel that contributed data to this project throughout their careers.

REFERENCES

Author notes