Relationship between Thermal Variance and Hatchery Estimates of Fertility in a Multispecies Assemblage of Anadromous Salmonids in Response to Managed Flows

Robert M. Sullivan ^1,*, John P. Hileman ²

1. California Department of Fish and Wildlife, Region 1, Wildlife/Lands Program, P.O. Box 1185 Weaverville, California 96093

2. California Department of Fish and Wildlife, Region 1, Fisheries Program, Trinity River Project, P.O. Box 1185, Weaverville, California 96093

* Correspondence: Robert M. Sullivan

Academic Editor: Petros Ganatsas

Special Issue: Environmental Monitoring and Protection of Inland River

Received: November 15, 2021 | Accepted: February 13, 2022 | Published: February 22, 2022

Adv Environ Eng Res 2022, Volume 3, Issue 1, doi:10.21926/aeer.2201009

Recommended citation: Sullivan RM, Hileman JP. Relationship between Thermal Variance and Hatchery Estimates of Fertility in a Multispecies Assemblage of Anadromous Salmonids in Response to Managed Flows. Adv Environ Eng Res 2022; 3(1): 009; doi:10.21926/aeer.2201009.

© 2022 by the authors. This is an open access article distributed under the conditions of the Creative Commons by Attribution License, which permits unrestricted use, distribution, and reproduction in any medium or format, provided the original work is correctly cited.

1. Introduction

Numerous laboratory and in-river field studies have demonstrated that thermal variance in water temperature is a critical factor affecting all aspects of adapted life history expressions in fish (Table S1). These studies demonstrate the potential for both direct and indirect effects of extreme fluctuations in water temperature on fertility, egg development, growth, timing of migration, timing of spawning [1]. Natural thermal regimes vary widely in complex and dynamic geomorphologically and topographically diverse river systems in both time and space. However, although natural environmental influences determine river flow and its effect on water temperature in unregulated rivers, these variables may be greatly altered by hydroelectric dam operations that redistribute water and regulate flow in highly managed riverine systems [2,3]. Anthropogenically-induced variability in hydrology and its impact on the thermal regimes of major rivers is a critical factor affecting the physiology, ecology, and life history strategies of resident salmonids [3,4,5,6,7,8]. The timing, variance, and magnitude of reservoir water releases have had dramatic effects on temperature regimes in watersheds and highly managed river systems [9,10,11]. Knowledge of the thermal regime of a river is critical to effective management of environmental flow, as well as understanding the potential impacts of dam-generated hydroelectric power on freshwater ecology and conservation of riverine fisheries resources; a topic recently reviewed and illustrated by Geist [3].

Water temperature models, with and without in-stream flow requirements suggest a clear relationship between river water temperatures and rate of flow [12]. This phenomenon directly influences the incubation temperature regime, which in turn influences the phenology and survival in embryos and alevins ([13]; Table S1). For example, Weinhemier et al. [14] showed that peak flows during incubation were strongly and negatively related to egg survival, and warmer temperatures were related to early emergence and timing of migration. Beckman et al. [15] identified genetic changes that affect timing of emergence under a common thermal regime, in which a genetic response and local adaptation to riverine thermal regimes may result from both mean and variance measures in water temperature. Similarly, highly modified natural flow regimes may cause life history strategies of fish assemblages, historically adapted to regional habitat templates, to be significantly altered downstream of dams in only a few decades owing to anthropogenic modifications in timing of stream flow [16,17,18,19]. Because salmonid populations have evolved in response to the temperature regime they historically used for migration and spawning, deviations from the normal pattern can affect survival [20]. Species with a narrow thermal ‘‘niche’’ are particularly vulnerable to alterations in water temperature regimes [21]; and taxa adapted to cold-water habitats may be especially vulnerable to the increases in temperature that commonly result from human influences [22,23]. Thus, it is critical that fisheries resource managers understand how variation in thermal physiology yields population-specific thermal performance [24].

Significantly, temporal variability in riverine water temperature has important biological and ecological consequences as a function of nonlinear relationships in species-specific physiological responses [25,26,27]. Yet adaptive management and conservation planning frequently only consider the effects of spatial heterogeneity in thermal variance. These efforts rarely reflect on the complexity of temporal heterogeneity [28,29] and its effect on the complex interactions between abiotic and biotic factors that relate time and resource use [30,31]. While spatial priorities generally derive from first-order temperature metrics (i.e., use of only mean values) traditional degree-day models are mostly insufficient to explain the effects of mean water temperature on the timing of life history strategies at regional scales [12,27]. Evidence continues to mount that biological and ecological reactions to thermally fluctuating environments will differ from expectations based on mean temperature alone [26]. This is because conventional analyses tend to oversimplify modeling by ignoring critical elements of thermal diversity by placing greater emphasis on simplistic lethal maxima underlying most water quality standards [32]. For example, in fisheries management use of landscape models to predict future distributions of species generally employ temperature means but not fluctuations in thermal variability [24,33,34,35]. Accentuating this complexity are alterations in the “natural” thermal variance originating from extrinsically managed hydrographs at seasonal, weekly, and daily scales as a function of reservoir operations, dam discharges, water diversions, and actions associated with in-river geomorphological restoration [3,36].

Knowledge of fish productivity in altered riverine systems is particularly important for management of anadromous salmonids given that native stocks have experienced significant declines due to regional environmental and anthropogenic factors [37,38,39,40,41]. Throughout the Central Valley of California, populations of anadromous Chinook Salmon and steelhead have experienced dramatic declines in the last century owing to construction of many water impoundments and diversions on many riverscapes. These alterations changed the historic flow and thermal regimes, resulting in thermally coupled deviations in fish development, growth, health, distribution, and survival [9]. In northern California, discharge of flows from Lewiston Dam into the upper Trinity River involves complex issues affected by political, economic, environmental, and biological factors. Depending upon the configuration of the hydrograph the thermal characteristics of the discharge can vary widely and must be considered relative to the potential effects on all aquatic biota [42]. In the upper-most cold-water reach of the Trinity River accessible to spawning salmonids there have been considerable and consistent fluctuations in the volume of flow since 2003. Consistent annual alterations of the “normalized” post-dam flow are in large measure a function of in-river managed geomorphic restoration actions [43]. Recent time series modeling suggests that subcomponents of the historical post-dam flow volume, keyed to geomorphologic restoration actions, significantly misrepresent the time series characteristics, seasonal-trend decomposition patterns, and forecast accuracy compared to the baseline pre-restoration time series hydrographic model [36]. These activities also have altered run-timing in anadromous hatchery- and natural-origin spawning salmonids Chinook Salmon, Coho Salmon, steelhead, and local migration patterns of non-anadromous Brown Trout and Klamath Small-scale Suckers relative to the historical post-dam pattern [44,45,46]. All of these species of salmonids spawn predominantly within the first 15 km of the mainstem below Lewiston Dam; and Chinook Salmon now increasingly spawn down-river and outside newly constructed mainstem restoration sites [47].

The purpose of our study focused on evaluating the relationship between thermal variance in river water and hatchery-measured estimates of female fertility as a function of anthropogenic altered annual flow. The question we pose is this. “What is the relationship between estimates of fertility in artificially propagated hatchery parental brood-stock and thermal variance in river water as a function of highly managed annual flow regimes?” To address this question, we: 1) describe and evaluate annual and seasonal variation in flow volume and water temperature as a function of continuous flow management since 2003; and 2) assess the pattern and relationship between flow-mediated water temperature and fertility within a multi-species assemblage of co-occurring anadromous salmonids. Understanding the relationship between flow management, temperature, and measures of fertility may help researchers improve their predictions of how temperature affects various aspects of species-specific phenology [7,48] in salmonid life history data within both a controlled hatchery environment and an in-river setting.

2. Materials and Methods

2.1 Study Area and Background on Managed Flows

Our study focused on the upper-most reach of the Trinity River in northwest California, USA (Figure 1). Trinity and Lewiston dams were constructed in the early 1960s. Trinity Dam creates Trinity Reservoir that stores up to 3,022 m³ of water. Lewiston Reservoir formed by Lewiston Dam is located 11.8 km downstream of Trinity Dam (river kilometer [rkm] 179.8). It functions as a re-regulating reservoir for flow into the Trinity River and as a source of water diverted to the Sacramento River Basin [49]. Lewiston Dam is the uppermost limit of anadromous fisheries on the Trinity River. Administratively, the upper Trinity River (mainstem) is the section of the river from the confluence of the North Fork Trinity River 63.1 rkm upstream to Lewiston Dam. The Trinity River Hatchery (henceforth called “hatchery”) is immediately below Lewiston Dam.

Figure 1 Map showing Trinity County (shaded black area in CA map), Trinity and Lewiston reservoirs and dams, Trinity River Hatchery; the upper-most limit of the primary traditional natural area spawning ~ 1.7 rkm downstream of the hatchery, including side-channel (dashed gold-colored line); and the Lewiston Water Quality gauge (LWS; https://earth.google.com/web). Black bars in inset map indicated where Trinity and Lewiston dams are located.

In 2000, the United States Department of Interior (USDI), Record of Decision (henceforth called the “ROD”) outlined a plan for restoration of the mainstem Trinity River and its fish and wildlife populations. The ROD was based 20 years of studies of the Trinity River and its fishery resources that culminated in the Trinity River Flow Evaluation Final Report [43,49,50], which recommended that high summer flows were necessary to maintain or restore salmon habitat and populations on the Trinity River. In 2000 the ROD directed that ~ 50% of the water would remain in the river (not be diverted to Central Valley). As a result of the ROD and the TRFES, the Trinity River Restoration Program [43] was formally established in 2002 to recommend how water was to be released for restoration of the mainstem Trinity River and its fisheries resources. The ROD defined five water year-types along with a minimum volume of water to be released into the Trinity River annually in accordance with the annual flow recommendations. The term “water year” defined by the United States Geological Survey is “the 12-month period 1 October (of any given year) through 30 September of the following year” (https://water.usgs.gov/nwc/explain_data.html).

The strategy for restoration of the mainstem Trinity River included: 1) flow management through manipulation of the annual hydrograph, 2) mechanical channel rehabilitation, 3) sediment management, 4) watershed restoration, 5) infrastructure improvements, 6) adaptive environmental assessment and monitoring, and 7) environmental compliance and mitigation [43]. These actions were intended to clean spawning gravels, build gravel bars, scour sand out of pools, provide adequate temperature and habitat conditions for fish and wildlife at different life stages, control riparian vegetation, and perform other ecological functions. From 2003 through 2019, total restoration releases have included flows for: 1) restoration, 2) late summer pulse flows, and 3) Tribal Ceremonial Boat Dances (Figure 2A and Table S2). Tribal flows occur only in odd years, just prior to any pulsed augmentation flows. Since 2003 they have amounted to ~ 0.65% of the total annual release into the Trinity River. Pulsed augmentation flows generally tier-off the trailing ends of ceremonial flows [44,45,46].

Figure 2 A) Representative examples of hydrographs of the three different flow-types (PreROD [2002], ROD [2005], ROD + Pulse [2015]). Tribal Ceremonial Boat Dance flows occur in odd years. Hydrographs associated with ROD and ROD + Pulse flows were typically implemented between Julian week JW13 and JW40. Shaded colored areas correspond to an overlay of the run-time presence of each species in the upper river for each flow-type example. B) Historical seasonal (month / JW) run-time distribution as determined by total hard-counts of Coho Salmon (1990 - 2018), spring- and fall-run Chinook Salmon (1994 - 2018), and steelhead (1990 -2018) obtained from Trinity River Hatchery.

2.2 Characteristics of Baseline and Managed Hydrographs

There are no comprehensive or systematic baseline data for pre-dam hydrographic or water temperature for the Trinity River below Lewiston Dam. Therefore, we equate “baseline” post-dam flows (henceforth called “Pre-ROD flows”) to flow releases that occurred after such data became available from 1994 to the present. Pre-ROD flows (1994 - 2018) are those that occurred after construction of the Trinity and Lewiston dams yet prior to implementation of annually managed hydrographs or pulsed flow augmentations, except on rare occasions due to an “emergency-of-dams” release because of abundant rain in combination with snow melt, as occurred in 1974. Pre-ROD flows generally lasted ~ 52.4 days from late April to late July (range 28.0 - 81.0 days) and averaged ~ 119.6 m³/second (range 62.3 - 192.3 m³/sec) of flow at the peak of the hydrograph (Figure 2A; [44,45]). In contrast, from 2003 to 2018 ~ 25.7 km³ of water was released from Lewiston Reservoir; of which 50.6% was released down the Trinity River and 49.4% was diverted to the Sacramento River. Of flow released down the Trinity River ~ 12.4 km³/year (x̄ = 0.77 km³) was allocated for in-river geomorphic restoration purposes (Table S2;[46]).

This restoration flow regime (henceforth called “ROD flows”) generally occurs from late April to August (Table S2). Timing, magnitude, shape, and duration of each ROD flow hydrograph varies annually depending upon the specific intent of the management action [43]. Average duration of each ROD flow was ~ 89.8 days from mid-April to early August (62.0 - 112.0 days) and averaged ~ 221.9 m³/second (124.9 - 328.6 m³/sec) of flow at the peak of the hydrograph [46]. Supplemental late summer pulse-flow augmentations have also occurred in 2003 and 2004, from 2012 to 2016, and in 2018. Implementation of these late summer pulsed flow augmentations were intended to cue upstream seasonal migration of Chinook Salmon out of the lower Klamath River to potentially reduce risk of density dependent epizootic of the ciliate parasite Ichthyophthirius multifiliis. Average duration of pulsed flows approximated 28.3 days from mid-August to late September (11.0 - 40.0 days) and averaged ~ 61.1 m³/second (35.3 - 97.0 m³/sec) of flow at the top of the hydrograph (Figure 2A; [46]). Importantly, we note that previously each pulsed flow event was also accompanied by a single ROD flow hydrograph (henceforth called ROD + Pulse flows) beginning in 2003. Thus, for each Pulse + ROD flow, the effects of each pulsed augmentation are not totally independent from the effects of its companion ROD flow.

2.3 Data Collection

Representative examples of each annual flow-type and distinguishing characteristics, and the seasonal run-time schedule specific to each species based on hatchery hard-counts of individual fish are illustrated in Figure 2A and Figure 2B. Except for analyses of annual variation (1994 - 2018) in each environmental parameter, all species-specific comparisons involving environmental and reproductive data were collected only during the run-time schedule that each taxon was historically in the upper river and not at other times (Figure 2B). In-river water temperature was obtained at the Department of Water Resources Lewiston Water Quality gauge (LWS; Figure 1), 1.7 rkm downstream from the Lewiston Dam and hatchery (rkm 178.2; 40°43’ 13N, 122°47” 28W; elevation 558 m). At the LWS gauge daily recordings of telemetered digital data were used to assess the seasonal extent of variability in the hydrology and thermal regime that characterized the upper Trinity River. Digital data included averaged daily: 1) maximum air temperature (ADAIR °C), 2) flow volume (ADFV m³/sec), and 3) water temperature (ADWT °C). Although punctuated up-and-down measurements using ADWT may approximate a bi-monthly framework for viewing this metric, it does not identify or illustrate extremes in daily water temperatures. Instead, the range (maximum - minimum) of extreme variability in daily water temperature (ADWTV °C) was used to address the need for such a metric [26,36].

Data obtained from the LWS gauge was used for several specific reasons. 1) There are no other sources of flow augmentation from the watershed or major tributaries between the LWS gauge and hatchery. So, this section of the river is not inadvertently influenced by conditions attributable to the local watershed. 2) Water temperature in the hatchery is monitored daily from the main intake throughout the facility to the river effluent, which historically is assumed to be within ~ 1°C (D. Muir, personal communication). Thus, the LWS gauge provides the most consistent and detailed long-term digital information reflecting water temperature used during each hatchery phase of fertilization, fertility estimation, egg production, juvenile grow-out, and release of hatchery-produced anadromous salmonids. 3) Virtually all hatchery-raised Chinook Salmon spawn within the first 3.2 km below the Lewiston Dam, which consistently contains the largest densities of hatchery-origin Chinook Salmon spawning redds [47]. 4) The LWS gauge is the closest gauge to the hatchery and provides the best location for measuring water temperature conditions for in-river spawning hatchery-origin Chinook Salmon and other anadromous species spawning in this segment of the reach. 5) The gauge provided data used in all National Environmental Protection Act (NEPA) assessments of flow augmentation and fluctuations in riverine average daily water temperature, specific to the upper Trinity River [35,51].

2.4 Timing and Measurement of Fertility at Hatchery

Our analysis of hatchery-collected data included all years for which hard-counts of hatchery and “wild” stocks of adult fish (> 32 cm fork length) were available for spring- and fall-run Chinook Salmon (1994 - 2018), Coho Salmon (1990 - 2018), and steelhead (1990 - 2018; Table 1). This process was conducted in accordance with a Julian week (JW) schedule, defined as one of seven consecutive-day-sets of 52 weekly periods in a calendar year that allowed inter-annual comparisons of identical weekly periods. Brood-stock were collected by fyke trap at the hatchery during annual egg-take operations, which generally last for six months out of the year (JW36 [Labor Day in early September] - JW11 [17 March]; Table 1). Spring-run Chinook Salmon egg-take starts in early September and goes into mid-October (JW36 - JW41). A break in spawning of hatchery Chinook Salmon (spawning break) occurs in JW42 to segregate races. Fall-run Chinook Salmon egg-take starts in late October and lasts to early December (JW43 - JW50). Coho Salmon egg-take starts in early November and lasts to mid-December (JW45 - JW51). Steelhead egg-take consists of 12 consecutive Tuesdays beginning at the end of the year (JW52 - JW11).

Table 1 Species-specific seasonal schedule of egg-take at the hatchery during annual Julian week (JW) artificial spawning cycles, including the number of available species-specific years sampled for reproductive data. Measurements of all reproductive variables were obtained at the time of egg-take at the hatchery during the run-timing that each species was in the river.

Measurements of reproductive variables were obtained for each species at the time of egg-take following a specific schedule (Table 1). Data were collected from fish that were ripe on the day of collection, including: 1) total number of eggs (ova) harvested at fertilization per spawn-day lot (TLEG) and 2) total number of eggs harvested that survived to the eyed ova stage or “final pick” (EAPIC). Dead eggs were removed daily, and each spawn-day lot was followed through to the final pick. Averaged daily percent fertility (ADPF) was calculated as: EAPIC / TLEG) * 100, which equates to percent of egg survival to the eyed ova stage for that spawn-day lot. Water was supplied to the hatchery from the base of Lewiston Dam. All developing ova were subject to ambient water temperatures that were within optimal species-specific thresholds. Traditional hatchery estimates of ADPF historically were used as an “index of fertility” derived from artificially propagated salmonids collected during the annual season when reproductively active adult cohorts of each species were in the upper Trinity River. Hatchery generated fertility estimates represent average spawn-day measurements based on the number of fish processed that spawn-day. The hatchery began supplementing the water in the spawning pans for fish eggs, with salt in the 2015 season, at a concentration of one cup per 0.24 liters of water. The salt solution was pre-mixed and added to each spawn-pan prior to sperm and eggs added. The intent was to reduce egg loss (i.e., improve fertility) in response to declining fertility in Chinook Salmon eggs brood-years 2012 to 2014 (personal communications L. Glenn [2014], and S. McCarn and S. Ballard [2015,2017]). We found no documents indicating that this procedure was modified prior to 2014; thus, our regression models only consider annual trends in fertility through 2014 for each species.

2.5 Statistical Analysis and Theoretical Distribution

All analyses were performed used R [52] and statistical significance was set at α < 0.05. Evaluation of univariate normality of temperature and flow volume measured at the LWS gauge was assessed using the adjusted Anderson-Darling (AD.adj) statistic [53]. This analysis suggested that all environmental variables were not normally distributed (ADFV: AD.adj = 1,480, p < 0.001, n = 8,401; ADWT: AD.adj = 8.7, p < 0.001, n = 8,055; ADWTV: AD.adj = 276.1, p < 0.001, n = 8,042). Hatchery-derived species-specific reproductive attributes also deviated significantly from normality (TLEG: AD.adj = 73.8, p < 0.001, n = 1,435; EAPIC: AD.adj = 116.1, p < 0.001, n = 1413; ADPF: AD.adj = 51.79, p < 0.001, n = 1,187). Thus, all statistical analyses of both environmental and reproductive attributes used non-parametric methods [54]. Kruskal-Wallis Chi-square rank sum tests (ꭓ²) evaluated each variable among each designated flow-type group. When the omnibus test was rejected post-hoc planned pairwise comparisons between each group were made using the Dunn test statistic (Z). P-values were adjusted using the Benjamini-Hochberg method [55]. The Spearman rank correlation coefficient (r_s) 2-tailed test calculated the strength and direction of the relationship between environmental and reproductive attributes [56]. Visual assessment of each environmental variable was conducted using two goodness-of-fit plots to determine the appropriate probability distribution (Figure S1; [57]). Q-Q plots and the Akaike information criterion goodness of fit statistic (AIC; [58]) suggested that a normal probability distribution more closely approximated ADWT. Whereas a Gamma probability distribution most closely approximated ADFV and ADWTV. Thus, the Gamma distribution was preferred in all subsequent regression analyses because it more correctly described the center and tails of the empirical distributions for each variable.

Semi-parametric Generalized Additive Models (GAM) were used in all regressions [59,60]. Response curves generated from each GAM illustrate the relationship between the fitted function and the response variable. Statistics reported from each GAM were: 1) GAM-F (approximate significance of smooth terms) and p-value, 95% confidence bands for spline lines [61], and adjusted regression coefficient (R² adj.) for the model. The regression coefficient explained the degree to which each predictor variable(s) (ADWT and ADWTV) explained the greatest amount of variation in the predicted variable (ADPF). The Spearman’s rank correlation coefficient was used to assess the strength and significance of each trend delineated by the smooth term [62]. Stationarity in each environmental variable (ADFV; ADWT; ADWTV) and estimate of fertility (ADPF) was evaluated in response to each species-specific flow-type hydrograph. Autocorrelation analyses of GAM-generated residuals were used to investigate the relationship between each time point relative to each previous time point in the distribution of lag values for each metric.

This process facilitated a test of the hypothesis that flow-mediated thermal variance in river water could potentially affect estimates of fertility in co-occurring species in a similar way. The Augmented Dickey-Fuller (ADF) statistic, Box-Ljung test (ꭓ²), and Breusch-Godfrey test (LM; [63]) were used to evaluate stationarity, and test for low- and high-order serial correlation across lags, respectively. Generic k-fold cross validation was used to estimate the strength and ability of our models to predict future observations [64]. Several different metrics were used to measure the quality of each model. 1) A t-test measured the overall significance of the model. 2) R² measured the strength of the linear relationship between the predictor variable(s) and the response variable; where an R² = 1.0 indicates a perfect linear relationship, while an R² = 0 indicates no linear relationship. 3) Root Mean Squared Error (RMSE) measured the average prediction error made by the model in predicting the value for a new observation. 4) Mean Absolute Error (MAE) represents the average absolute difference between the true value of an observation and the value predicted by the model; it is generally less sensitive to outliers compared to RMSE. For both RMSE and MAE lower values indicate a better fit to the model.

2.6 Ethics Statement

Research presented here did not involve human, animal, or plant subjects, cell lines or clinical trials. All analyses were conducted with data previously collected at the hatchery or generated from physical parameters collected from the riverine environment.

3. Results

3.1 Annual and Seasonal Trends in the Historical Time Series of Environmental Variables

Significant overall differences in the historical annual trends of each flow-type hydrograph were observed for each environmental variable measured in-river at the LWS gauge (Table 2). For ADFV the baseline PreROD flow-type exhibited the greatest amount of variation relative to other flow-types as measured by its standard deviation, which reflected a more climatically “natural” post-dam hydrograph compared to regularly manipulated managed hydrographs. For ADWT there was a distinct transition in the increased level of variance from PreROD and ROD to ROD + Pulse flows. All GAM regressions suggest that smooth terms for each environmental variable were non-linear and significant for each annual response curve, but the amount of variation explained by all regression coefficients was not robust ranging from R² adj. = 0.3% (ADFV) to R² adj. = 17.8% (ADWTV; Figure 3A, 3B, and 3C). Exaggerated variance in ADFV and ADWT was also evident in higher volumes of flow in years of ROD and ROD + Pulse flows. This pattern was particularly evident in extreme values of thermal variance (ADWTV). A clear distinction between ROD and ROD + Pulse flows on an annual basis was observed in 2017 for both ADWT and ADWTV (Figure 3B and 3C).

Table 2 Kruskal-Wallis Chi-square rank sum tests (ꭓ²) of each environmental variable measured at the Lewiston Water Quality (LWS) gauge for each managed flow-type hydrograph (PreROD, ROD, ROD + Pulse), followed by planned post-hoc Dunn tests (Z) of all pairwise comparisons. Measurements of all flow and water temperature data were obtained at the Lewiston Water Quality gauge during the run-timing that each species was in the river. Included were the arithmetic mean, lower (LCL) and upper (UCL) 95% confidence limits, and the standard deviation (SD) for averaged daily: 1) flow volume (ADFV m³/sec), 2) water temperature (ADWT °C), and 3) extreme variability in water temperature (ADWTV °C). Probability levels were calculated at: < 0.05 = *, < 0.01 = **, < 0.001 = ***.

Figure 3 Composite boxplots of annual variation in average daily: A) flow volume (ADFV m³/sec), B) water temperature (ADWT °C), and C) extreme variability in water temperature (ADWTV °C). Seasonal Julian week (JW) fluctuations are illustrated in: D) ADFV, E) ADWT, and F) ADWTV. Shaded smooths (colored yellow) indicate 95% point-wise standard error for each curve surrounding each fitted Generalized Additive Model (GAM) function (centered black line) in A - B and centered red line in D - F). GAM regression statistics and Spearman rank correlation coefficients correspond to daily measurements of each variable from 1994 through 2019.

Here ROD flows resulted in lower values for both water temperatures variables (cooler and less variable) in conjunction with higher volumes of flow compared to ROD + Pulse flow years (2005 - 2011, 2018) immediately preceding and after pulsed augmentation flows in 2018. This pattern is consistent with the seven consecutive ROD flows from 2005 to 2011. Ranked correlation analyses found significant positive relationships between year and both water temperature variables but the correlation between year and ADFV was negative. As expected, the correlation between ADWT and ADWTV was significant and positive (r_s = 0.52, p < 0.001, n = 8,042), but both water temperature variables were negatively correlated with ADFV (r_s = -0.09, p = 0.001, n = 8,012 and r_s = -0.27, p < 0.001, n = 7,999, respectively) and the strengths of the correlations were weak. Overall, patterns depicted by the regression and correlation analyses suggested that the historical water temperature in the upper-most reach below the hatchery tended to increase slightly on an annual basis, while the corresponding volume of flow tended to decline for the years sampled.

Similarly, the historical trends expressed in seasonal Julian week fluctuations in each environmental metric also exhibited non-linear patterns (Figure 3D, 3E, and 3F). Again, although all smooth terms were significant for each GAM response curve the relationship between the predictor (Julian week) and each environmental response variable was weak; with R² adj. ranging from 58.5% (ADWT) to 16.2% (ADWTV). There was a significant positive relationship between Julian week and ADWT, but this relationship was negative for both ADFV and ADWTV. Although a seasonal component in the thermal regime extending from winter to fall was evident in ADWT, seasonal hydrographs describing ADFV were augmented by managed hydrographs. This pattern was reflected by extreme thermal variation in water temperature (ADWTV) as mirrored by the correlation with Julian week. Extrinsic to the thermal regime of flow volume we found that along this section of the upper mainstem the relationship between maximum air temperature (ADAIR) and both water temperature variables was positive and significant (ADWT: r_s = 0.672, p < 0.001, n = 5,664 and ADWTV: r_s = 0.48, p < 0.001, n = 5,664). The weak relationship between ADAIR and extreme thermal variance in water temperature, despite large sample size, may be a function of the nearness to the outlet at the base of Lewiston Dam in combination with the insulating effect of steep topography and narrowness of the adjacent canyon walls along this section of the river.

The extent of seasonal differentiation and dispersion among flow-types for each environmental variable was further evident when Julian week data were evaluated by flow-type (Figure 4). This process revealed seasonal irregularity in more traditionally cyclical fluctuations of each environmental variable for PreROD flows compared to: 1) high-volume ROD flow hydrographs employed systematically and continuously for in-river geomorphologic restoration (JW18 - JW22); and 2) ROD + Pulse flows intended for both in-river geomorphological work and secondarily to prompt upstream seasonal migration of Chinook Salmon out of the lower Klamath River to reduce the potential for density dependent disease transmission. Pulsed augmentations in conjunction with their companion ROD flow, and in odd years ceremonial flows (JW36 - JW40), accounted for more Julian weeks with extreme levels of thermal variance than either PreROD or individual ROD flows (Figure 4B, 4E, and 4H). This response was evident also by the saddle-shaped distribution of ADWTV on both sides of JW20 (Figure 4E[b], 4H[b], and 4I[b]), which was indicative of reduced thermal variability at the top of each annual hydrograph in volume of flow (Figure 4B[b] relative to 4H[b]). The greatest amount of variability in water temperature was associated with both the ascending and descending slopes of peak discharge during ROD and ROD + Pulse flows.

Figure 4 Flow-type boxplots of seasonal Julian week (JW) variation in average daily: flow volume (ADFV; A - C), water temperature (ADWT; D - F), and extreme variability in water temperature (ADWTV; G - I). Lower case lettering reference subsections of each seasonal flow-type described in the text. Shaded smooth (colored yellow) indicate 95% point-wise standard error for each curve surrounding each fitted GAM function (centered black line); JW arranged sequentially beginning with JW01 (1 - 7 January) to JW52 (24 - 31 December). GAM regression statistics correspond to daily measurements of each variable from 1994 through 2019.

Importantly, the extreme flow-mediated water temperature fluctuations described above for each environmental variable resulting from managed hydrographs post-2002 were not characteristic of baseline PreROD flows. Similarly: 1) measurements obtained in-river at the LWS gauge indicated significant overall differences between flow-types for each environmental variable when evaluated in accordance with the historical run-time schedule of each species (Table 3); 2) planned post-hoc tests showed that 62.5% (n = 48) of all pairwise comparisons were significantly different among flow-types for each species-specific environmental variable; and 3) the standard deviations for all environmental variables were consistently higher for ROD + Pulse flows and lowest for baseline Pre-ROD flows in all species comparisons.

Table 3 Kruskal-Wallis Chi-square rank sum tests (ꭓ²) of environmental variables measured at the Lewiston Water Quality gauge for each managed flow-type hydrograph (PreROD, ROD, ROD + Pulse) for each species of salmonid, followed by planned post-hoc Dunn test statistic (Z) of all pairwise comparisons. Measurements of all flow and water temperature data were obtained at the Lewiston Water Quality gauge during the run-timing that each species was in the river. P-values adjusted using Benjamini-Hochberg method. Included are the arithmetic mean, lower (LCL) and upper (UCL) 95% confidence limits, and standard deviation (SD) for averaged daily: 1) flow volume (ADFV m³/sec), 2) water temperature (ADWT °C), and 3) extreme variability in water temperature (ADWTV °C). Probability levels: < 0.05 = *, < 0.01 = **, < 0.001 = ***.

3.2 Variation in Reproductive Attributes of Hatchery Parental Brood-Stock

Annual hatchery production goals are set by the United States Bureau of Reclamation. To meet these goals the hatchery egg-take (TLEG) for spring-run Chinook Salmon was set at 3 million eggs, fall-run Chinook Salmon at 6 million eggs, Coho Salmon at 1.2 million eggs (reduced to 720,000 in 1990), and 2 million eggs for steelhead (reduced to 1,200,000 in 1990; Figure 5A - 5D). However, these goals: 1) have not been met since 2014 for spring-run Chinook Salmon; 2) only attained in 2018 since 2012 for Fall-run Chinook Salmon, 3) have not been met since 2015 for Coho salmon, and 4) have not been met since 2016 for steelhead. Alternatively, if EAPIC is used as an index the annual egg take of fertilized viable eggs for each species was much lower than for TLEG (Figure 5E - 5H). In contrast, TLEG harvested for Spring-run and fall-run Chinook and Coho salmon suggested negative trends for the same period (r_s = -0.29, p = 0.001, n = 290; r_s = -0.15, p = 0.004, n = 392 and r_s = -0.20, p = 0.001, n = 395; respectively). Whereas the pattern in steelhead suggested a significant positive trend (r_s = 0.47, p = 0.001, n = 342). Except for steelhead averaged spawn-day trends in TLEG harvested for the years sampled reflected patterns consistent with sharp declines in stocks of all taxa since 2003 [45].

Figure 5 Bar graphs summarizing the total number of eggs harvested (TLEG; colored blue) and the total number of eggs left after the final pick (EAPIC colored red) by year from adult female salmonids spawned at the Trinity River Hatchery by species: spring-run Chinook Salmon and fall-run Chinook Salmon (1994 - 2018), Coho Salmon (1990 - 2018), and steelhead (1990 - 2018). Years with maximum numbers of eggs harvested are indicated for both variables and dashed horizontal lines (colored black) correspond to species-specific egg-take goals discussed in the text.

For all taxa except Coho Salmon there were significant overall differences in EAPIC (Table 4). Similarly, there were significant species-specific differences in 75.0% (n = 12) of all planned in post-hoc pairwise comparisons between flow-types. Except for Coho Salmon, the largest discrepancies between TLEG and EAPIC were observed in managed flow-types (Figure 6A - 6D); while the smallest discrepancies were found in baseline PreROD flows. These trends also were reflected by consistently higher variances (as measured by SD) in each species-specific flow-type for ROD and ROD + Pulse flow hydrographs compared to baseline PreROD flows. Accordingly, for Chinook Salmon and steelhead, irrespective of run-size and oversampling of eggs by the hatchery to meet mandated quotas, our results suggest that in years of managed hydrographs (to facilitate in-river restoration actions) there were significantly fewer viable eggs obtained at the hatchery compared to years of baseline Pre-ROD flows. Currently these subtle differences in estimates of fertility relative to flow-type are not being monitored.

Table 4 Kruskal-Wallis Chi-square rank sum tests (ꭓ²) of the differences between total number of eggs collected and total number of eggs left after pick measured at the Trinity River Hatchery between flow-type hydrographs (PreROD, ROD, ROD + Pulse) for each species of salmonid, followed by planned post-hoc Dunn test statistic (Z) of all pairwise comparisons. Included are the arithmetic mean, lower (LCL) and upper (UCL) 95% confidence limits, and standard deviation (SD). Probability levels: < 0.05 = *, < 0.01 = **, < 0.001 = ***.

Figure 6 Polygons of species-specific flow-types illustrating the differences between the total number of eggs harvested (TLEG) and the total number of eggs left after the final pick (EAPIC) to remove dead eggs (TLEG - EAPIC).

3.3 Autocorrelation in Species-Specific Annual Trends of Environmental Variables and Estimates Fertility

Autocorrelation analyses of annual variation in both flow and water temperature variables and estimates of fertility by use of GAM residuals suggested virtual unanimity in stationarity. All lags were contained within approximate 95% confidence limits for each variable in all species comparisons. All ADF-statistics were significant indicating stationarity in each comparison (Figure S2 and Figure S3). Both the Box-Ljung and Breusch-Godfrey tests indicated a lack of low- or high-order serial correlation among GAM residuals for all comparisons. Additionally, significant overall differences were observed in ADPF for all species-specific flow-types as indicated by Kruskal-Wallis Chi-square rank sum tests; and 83.3% (n = 12) all post-hoc planed pairwise tests suggested significant differences between flow-types for each taxon (Table 5). Again, results indicated that for all species except fall-run Chinook Salmon, the greatest amount of dispersion (SD) in estimates of fertility were found in managed flows, especially in ROD + Pulse flow-years. Prior to the change in the methodology for estimating in fertility beginning in 2015, all salmonids artificially propagated at the hatchery experienced a significant annual decline in fertility through 2014 (Figure 7A - 7D). Post-2014 estimates of fertility in all species except steelhead appeared to suggest an initial increase but in subsequent years appeared to decline through 2018. However, there are too few data to make a definitive statement about future trends in post-2014 estimates as methods for measuring fertility continue to be altered annually, which unfortunately prevents the possibility of developing a historical time series model that includes post-2014 estimates.

Table 5 Kruskal-Wallis Chi-square rank sum tests (ꭓ²) of the differences between species in averaged daily percent fertility (ADPF) estimated at the hatchery between flow-type hydrographs (PreROD, ROD, ROD + Pulse), followed by planned post-hoc Dunn test statistic (Z) of all pairwise comparisons. Included are the arithmetic mean, lower (LCL) and upper (UCL) 95% confidence limits, and standard deviation (SD). Probability levels: < 0.05 = *, < 0.01 = **, < 0.001 = ***.

Figure 7 Generalized Additive Model (GAM) regression plots of annual fluctuations in averaged daily percent fertility (ADPF) association with species-specific flow-types. Shaded smooth (colored gray) areas indicate 95% point-wise standard error distribution for each curve surrounding each fitted GAM function (centered black line). Annual hydrograph designations are indicated for baseline PreROD, ROD, and ROD + Pulse flows. GAM regression statistics and Spearman ranked correlation coefficients are for years 1994 to 2018 in both runs of Chinook Salmon and years 1990 to 2018 in each run of Coho Salmon and steelhead.

3.4 Fertility as a Function of Water Temperature and Flow Volume

Results of GAM regression showed a significant negative relationships between ADPF and both individual and combined water temperature variables (Table 6), and between flow-types for species-specific water temperatures (Table 3) measured in-river at the LWS gage for Chinook and Coho salmon, but not steelhead (Figure 8 and Figure 9). The negative relationship between ADPF and ADWT and between ADPF and ADWTV implies that as the temperature and thermal variability of river water increases estimates of fertility tends to decrease during species-specific run-time schedules. For example, in all polygon plots keyed to baseline PreROD flow, ADPF generally plotted above the regression line indicating relatively higher estimates of fertility in all species for this flow-type. Conversely, samples that identified managed flow-type hydrographs tended to plot below the regression line indicating relatively lower estimates of fertility for each species of salmonid. This was the situation whether the smoothed regression line was significant or not (ADWT for Coho Salmon and steelhead; Table 6). In all taxa this down-turn in the smoothed regression line for fertility typically occurred between ~ 7.5°C and ~ 9.0°C in ADWT in all species. The comparatively low end of the regression line in steelhead (~ 5.0°C) marked the beginning of the run-time spawning schedule in early spring in managed flow-type years (Figure 8D; Figure 2). A summary of species-specific cross-validation R² performance errors suggested a close match with results of the GAM (R²) regression models used to predict estimates of fertility (Table 6). Except for ADWT (R² = 35.8%) and ADWT + ADWTV (R² = 34.8%) in spring-run Chinook Salmon the error variances explained were relatively low in each model for each species; whereas error measures for both MAE and RMSE were relatively high and similar across all cross-comparisons. However, the regression coefficient for both the GAM and cross-validation models decreased in parallel going from spring-run Chinook Salmon to steelhead, which corresponds to the timing at which each species is in the upper Trinity River and the sequential onset of annually managed flows.

Table 6 Summary of species-specific Generalized Additive Modeling (GAM) and cross-validation model performance error measures for water temperature metrics (independent or predictor variables) versus average daily percent fertility (ADPF; dependent or response variable). ADWT = averaged daily water temperature; ADWT averaged daily extreme variability in water temperature; R2 = R-squared; RMSE = root mean-square error; MAE = mean absolute error. Probability levels: < 0.05 = *, < 0.01 = **, < 0.001 = ***.

Figure 8 Polygons of species-specific fertility (ADPF) and Generalized Additive Model (GAM) regression plots showing the relationship between ADPF and average daily water temperature (ADWT) by flow-type. Shaded smooth (gray colored) areas indicate 95% point-wise standard error distribution for each curve surrounding each fitted GAM function (centered black line). Vertical dashed black lines indicate water temperature ranges between ~ 7°C and ~8°C for each comparison at which a measured decline in fertility first became noticeable.

Figure 9 Polygons of species-specific fertility (ADPF) and Generalized Additive Model (GAM) regression plots showing the relationship between ADPF and extreme variability in average daily water temperature (ADWTV) by flow-type. Shaded smooth (gray colored) areas indicate 95% point-wise standard error distribution for each curve surrounding each fitted GAM function (centered black line).

4. Discussion

4.1 Implications of Flow-Regulated Water Temperature

Our study consists of several key findings. First, we found highly significant overall differences between flow-types in averaged daily flow volume, water temperature, and extreme variability in water temperature when evaluated in reference to the historical run-time schedule of each taxon. Second, we found significant differences among species-specific flow-types for each environmental variable. Third, we identified significant declines in the averaged daily number of total eggs harvested in all species except steelhead. Fourth, we provide evidence of significant overall differences in estimates of fertility for all species-specific flow-types, in which 83.3% of post-hoc planed pairwise comparisons revealed significant differences between flow-types for each taxon. Fifth, we observed significant annual declines in the historical estimates of fertility in all species through 2014 prior to a change in methodology for estimating fertility at the hatchery. That all species of spawned hatchery salmonids exhibited concordant declining trends in fertility suggests that taxa with divergent life history traits were responding to environmental conditions in parallel. When this measure of annual “fitness” drops in one species it also declined in all other sympatric taxa, whether in response to extrinsic events in the upper Trinity River or elsewhere. These negative species-specific annual trends in fertility were consistent with declining stocks observed in all species of anadromous salmonids that cohabitate the upper Trinity River [45]. Sixth, we documented significant negative correlations between species-specific estimates of fertility and both measures of water temperature in all taxa except steelhead. This may be that steelhead initiates spawning in late winter when the thermal regime is at its coldest before transitioning to warmer temperatures through the early fall spawning season. We suggest that a positive regression between fertility with increasing water temperature in steelhead is an expression of a diverse life history portfolio in which not all individuals are “genetically pre-programed to die” after spawning (semelparous) unlike salmon [65,66]. This evolved adaptation in steelhead may contribute to the species’ resiliency and persistence despite anthropogenic effects and environmental variability potentially affecting growth and survival [67,68,69,70,71].

Our study focused on evaluating the statistical relationship between anthropogenically altered flow on hatchery-measured fertility as a proxy for use in modeling reproductive success in an assemblage of endemic salmonids. This process is necessary in identifying annual trends in reproductive performance among a multi-species assemblage of fish with different life history adaptations. We did not test for a “causal” interaction between flow-mediated water temperature and in-river estimates of fertility in spawning of hatchery- or natural-origin fish below Lewiston Dam, because these data are unknown. Similarly, the potential for flow-mediated thermal variance and its effects on fertility also remains unknown for fish spawning further downstream and in proximity to newly constructed restoration sites since 2003. However, it is reasonable to hypothesize that the cumulative effect of flow-mediated thermal variance on estimates of fertility in salmonids spawning in-river would be much more dramatic than what we have documented for hatchery fish. This is because of the greater thermal variance owing to fluctuations in flow volume, geomorphology, topography, and extrinsic climatic conditions found along this river system compared to more “controlled” and spatially confined hatchery conditions. Use of estimates of fertility in stocks of artificially spawned hatchery fish would appear to be a practicable baseline index for: 1) monitoring species-specific annual productivity and individual performance; and 2) providing a measure of relative “fitness” within a controlled hatchery environment for same-time contrasts with in-river hatchery and natural-area spawning fish.

We acknowledge that correlation and regression analyses do not necessarily equate to causation. Nevertheless, our results suggest that extreme variation in water temperature, as a function of highly managed flow regimes, has the potential to negatively affect fertility in several anadromous taxa with adaptively divergent life history adaptations. We suggest that concordance patterns among sympatric anadromous species are useful in providing initial insight and confidence into understanding the variance relationships in these metrics [72]. As with all correlative analyses there are noteworthy caveats: 1) fertility estimates used in our analysis were from hatchery escapement fish only and egg-take was determined by a process that was 100% artificial; and 2) although trends in fertility addressed herein may affect “wild” populations, similar empirical estimates of fertility in hatchery- and natural-origin spawning salmonids do not currently exist for the upper Trinity River. This is because in-river spawning has historically been monitored solely through indexing by annual surveys of carcass counts.

We note that cross-validation errors trended high across all species-specific comparisons, suggesting that our GAM models may be somewhat inefficient in their applicability to use future “out-of-sample testing” or unknown data [72,73]. However, regression coefficients for both GAM and cross-validation models decreased in parallel going from spring-run Chinook Salmon to steelhead. Moreover, species-specific life-history strategies corresponded to the timing at which each taxon is in the upper Trinity River and the onset of annually managed flows. We hypothesize that model relevance and efficiency would be improved if in-river estimates of fertility and thermal variance were used during the run-time spawning schedules in which reproductively “ripe” fish were actually in the upper river. Importantly, both GAM regression and cross-validation models identified declines in estimates of fertility in relationship to water temperature going from spring-run to fall-run Chinook Salmon, Coho Salmon, and to steelhead. This pattern is plausible given that: 1) spring-run Chinook Salmon are more exposed to managed flow regimes on an annual basis than any other anadromous species because their presence in the upper river coincides with the onset of ROD and ROD + Pulse hydrographs, which are function of annual long-term geomorphological restoration actions in combination with pulsed flow events; and 2) the sequence and timing of spawning is such that fall-run Chinook Salmon, Coho Salmon and steelhead occur in ordered succession progressively later in the year when annual ambient water temperatures naturally drop.

4.2 Flow, Temperature, and Other Co-Variates

Our study sought to identify a regional explanation for declining fertility in stocks of co-occurring salmonids in the upper Trinity River. As such, we focused on the relationship between water temperature and anthropogenic altered volumes of flow. Yet other co-variates not yet studied in both salt and fresh water may also affect annual and seasonal patterns of relative abundance, timing of migration, and estimates of fertility in spawning salmonids. Some of these co-variants include genetics, nutrition in the marine life history stage, and water quality in the freshwater upstream migration life history stage. A current narrative is that returning adult Chinook in the Trinity River, as well as other stocks along the west coast, exhibit clinical signs of thiamine deficiency in the last two spawn cohorts [74,75]. Similarly, ocean conditions are frequently cited as the primary cause for decreasing stocks of anadromous salmonids in the Pacific Northwest as a function of food availability in both adult and juvenile fish in association with regional drought. However, yet to be fully documented are the linkages and the extent to which regional drought, thiamine deficiency, degraded watersheds, or conditions in the ocean individually or cumulatively pose direct threats to inland fisheries [38,40,76,77,78]. To be useful, information on these and other potentially pertinent co-variates must be known, available, reliably sampled, and consistently gathered throughout a defined sampling period. Such requirements are critical in time series modeling, model evaluation, and forecasting. Yet these data are not currently available for the Trinity River and most all other riverine systems. What we do known is that anthropogenic flow manipulation does influence the behavior of pre-spawn fish during migration if fish migrate too fast or too early [45,46]. For Chinook Salmon, the effects of flow on fish also impose severe logistical constraints for management of hatchery operations and compliance with established mandates, which potentially affect estimates of fertility by narrowed ripening curves and delaying gametic maturation. What is lacking is a comprehensive assessment of the effects of managed flows on the condition of pre-spawn fish in route to the hatchery and in the traditional spawning grounds downstream during species-specific run-time scheduling.

We also acknowledge that the potential effects of annually managed flows on thermal variance of the mainstem Trinity River occurs predominantly before spawning, except in early spring-run and early arriving fall-run Chinook Salmon [45]. This does not mean, however, that subsequent development to juveniles is unaffected by managed flows if extreme thermal variability is a potential factor affecting in-river spawning [Table S1]. For example, extreme fluctuations in water temperature, as a function of managed hydrology, coincide with the in-river run timing and the duration of life history events for all four taxa (i.e., ontogeny of fry emergence, juvenile rearing, and juvenile and smolt out-migration) in the Trinity River. These events are diagrammatically outlined in the life history schedule for each species in the Trinity River Flow Evaluation Final Report, which was written prior to initiation restoration actions [49]. Yet evaluation of the relationships between flow-mediated water temperature and the potential effects on estimates of fertility, behavior, or other life history attributes have not been the focus of in-river: 1) spawning and production assessments in relation to potential flow-fish-effects, 2) geomorphological restoration actions relative to baseline “control” conditions, or 3) investigations into other co-variates that also could affect fitness.

Pre-ROD flows pre-date (~ 60 years) post-dam construction and restoration. As such, pre-ROD flows are perhaps more similar to what endemic species of salmonids evolved with because high flow events more closely tracked weather conditions and pre-dam hydrographs. Collectively, our results suggest that these differentiated species are tracking thermal conditions in the upper Trinity River in parallel. This pattern may in part be a function of managed hydrologic events post-2003, which deviate dramatically from variability associated with the new “post-dam norm” [36]. Importantly, pre-restoration flows and thermal regimes that have characterized seasonal hydrographs in concert with regionally unimpeded tributary flow events automatically factor in annual variability in climate change. Given the continuation of emerging evidence of direct and indirect effects of thermal variance of river water on salmon fertility, growth, development, and susceptibility of salmonid eggs to in-river survival and altered emergence timing [25,77,79], it seems prudent that resource managers prioritize consideration of the potential effects of fluctuating thermal variance associated with anthropocentric flow manipulation in highly regulated rivers. This approach is particularly critical in relation to: 1) the distribution, abundance and habitat attributes of salmonid redds, 2) successful egg incubation and follow-on development to fry emergence, 3) outmigrant survival, and 4) quality of newly created juvenile rearing habitat as a function of altered floodplain geomorphology. The latter of which may adversely impact species-specific scheduling of regionally adapted behavioral life history strategies by resident taxa. Nevertheless, there are also other more proximate factors that influence the biology of resident salmonids at the regional and local levels. Increasingly, the lack of knowledge resulting from integrated, well-designed, and long-term restoration programs that have weak or questionable validity, is in large measure the basis of considerable criticism by the scientific community that river restoration actions are poorly designed and fail too often [80,81,82,83,84].

4.3 Management Considerations

Scientific and policy communities acknowledge the need for maintaining or restoring natural flow variability to sustain the ecological health of fluvial ecosystems given the growing evidence of the potential negative consequences of altered flow regimes on these systems and the fisheries they support [26]. Understanding the potential effects of altered flow regimes on fisheries productivity and the linkage to extreme fluctuations in the thermal regime facilitates avoidance of uncertainties when assessing flow-fish productivity to enable adaptive management [19,24]. Because alterations in the variance of riverine thermal regimes occur in response to anthropogenic events, these issues may inevitably become more problematic as climate change predictably progresses in-step with the demand for water [10,25,85,86,87]. In the mainstem Trinity River numerous reports have addressed fluctuations in water temperature, but only in the context of average daily water temperature [35,52,88,89]. This was the metric used to determine potential direct, indirect, and cumulative effects to the affected environment associated with late summer pulse flow that were supplemental to the ROD flow from Lewiston Dam in 2015 and 2016 [90]. Yet there has been no attempt to equate, measure, or monitor the potential effects of thermal variance on estimates of in-river fertility or reproductive output in anadromous salmonids associated with altered flow regimes initiated in 2003. Such an effort is crucial as late summer pulse flows temporally occur during the upstream migratory phase of the lifecycle of adult salmonids. Thus, contrary to the statement by Zedonis [91], the effects of augmented flow management on Trinity River stage and thermal variance are not well known, particularly as relates to a post-dam construction baseline time series model [36], potential effects on the ecology of the river system, and the life history adaptations of co-occurring anadromous salmonids.

Given the results and limitations of our study, we suggest that a more refined and long-term research focus is needed. Such a re-examination would allow simultaneous annual comparisons of fertility rates and egg-to-outmigrant survival at both the hatchery and in-river cohorts of the spawning population. The potential for thermal variance in river water, as a function of flow management, to effect biological and ecological processes suggests that more emphasis should be placed on collection, analysis, monitoring, and archiving variability in long-term flow and water temperature data. As noted by Zillig et al. [24] “Simple static temperature criteria can be improved by incorporating local data on salmonid fundamental physiology and on ecological conditions to produce population-specific thermal management strategies.” We believe this suggestion also applies in a coordinated way to both hatchery operations and simultaneously down river at multiple locations within traditional natural area spawning sites. Examination of the effects of thermal variance, seasonality, and heterogeneity at multiple time scales seems more appropriate than simply applying a mean value to monthly, weekly, or daily variation in water temperature and flow volume. A study design of this kind has been demonstrated successfully for Chinook Salmon within the Chiwawa River, Washington [92]. This recommendation strongly relates to valuation of in-river fertility, production in salmonid populations, and the need to address declining stocks of salmonids in conjunction with hatchery operations and restoration of the mainstem Trinity River [86].

As such, we offer the following recommendations. 1) Compile historical data that relate ecological responses to variation in flow and water temperature consistent with seasonal run-time schedules of focal species at the hatchery and simultaneously in-river to compare with measures of species-specific fitness and habitat quality [92,93,94]. 2) Initiate Science-based experiments to determine empirical responses by salmonids in response to timing and volume of flow, and thermal variance. 3) Use the un-inhibited headwaters above Trinity Reservoir and major tributaries of the Trinity River to model and gauge ecological response to flow. 4) Modify existing spawning surveys to evaluate redd and carcass data for specific geo-rectified spawning criteria [95,96] in order to better understand the underlying mechanisms of ecological response to flow in both the mainstem and major tributaries of the Trinity River. This approach has been used in assessing salmonid habitat suitability below Lewiston Dam [97] and juvenile rearing habitat at channel rehabilitation sites along the Trinity River [98]. 5) Monitor and model variation in the thermal regime at the hatchery from egg-take to juvenile outmigration in conjunction with same day collection of air and water temperatures, and flow data within traditional spawning grounds immediately below the inflow into the Trinity River from Lewiston Dam, and selected sites downstream within the mainstem corridor. 6) Annually randomly sample and monitor via regression analyses the relationship between fork-length (i.e., size distribution) and reproductive metrics (number of eggs, size, estimates of fertility) in representative subsamples of individual/age cohort species-specific hatchery spawners. 7) Evaluate habitat quality at the location of each carcass and redd occupied by live un-marked females that do not fit the category of “hatchery spawning” fish based on proximity to the hatchery [47]; and assess the reproductive condition of female whole carcasses that link biological data to habitat data gathered at each site.

Acknowledgments

We thank several anonymous reviewers for their helpful editorial comments, suggestions, and numerous probing questions, which in several instances will require future research despite our best efforts to address them.

Author Contributions

Robert Sullivan developed the concept, analyzed and modeled the data, developed the graphics, and wrote the manuscript. John Hileman provided the initial raw data, knowledge and discussion about hatchery operations and fish biology in general, and edited portions of the manuscript.

Funding

Funding was not provided for this research effort.

Competing Interests

The authors declare that no competing interests exist.

Additional Materials

The following additional materials are uploaded at the page of this paper.

1. Table S1: Examples of phenology, life stage, and life history characteristics affected by thermal variance in water temperature or flow volume referenced to citation and outcome.

2. Table S2: Water-year volume and flow-type summaries for managed hydrographs used for in-river restoration (ROD and ROD + Pulse flows) for Trinity River releases and diversions to the Central Valley of California since 2003. Modified from the Trinity River Restoration (TRRP 2021) webpage, which provides explanations as to how water volume metrics were calculated (https://www.trrp.net/restoration/flows/summary); ac-ft = acre feet; cfs = cubic feet per second in accordance with the webpage.

3. Figure S1: Goodness-of-fit plots for normal and Gamma variable distributions fitted to measures of averaged daily: A) flow volume (ADFV m3/sec). B) water temperature (ADWT °C), and C) extreme variability in water temperature (ADWTV °C). Plots include theoretical density (colored blue) and quantile-quantile (Q-Q) plots (colored red).

4. Figure S2: Plots showing time series autocorrelation analyses of environmental variables measured at the Lewiston Water Quality gauge for each species during the time that they were in the upper river. Analyses used residuals derived from Generalized Linear Models (GAM): A - C) spring-run Chinook Salmon (SRC), D - F) fall-run Chinook Salmon (FRC), G - I) Coho Salmon (COS), and J - L) steelhead (STL). Environmental variables were average daily: flow volume (ADFV m3/sec), water temperature (ADWT °C), and extreme variability in water temperature (ADWTV °C). Augmented Dickey-Fuller test = ADF, Box-Ljung test = ꭓ2, and Breusch-Godfrey test = LM. Horizontal axis shows different lag values, vertical axis shows the autocorrelation at each lag value, and 95% confidence limits (bands) for the ACF function correlogram (blue dashed horizontal lines).

5. Figure S3: Plots of the results of time series autocorrelation (ACF) analyses of averaged annual percent fertility (ADPF) obtained from fish spawned at the hatchery at the time that they were in-river. Analyses use residuals derived from Generalized Linear Models (GAM) for each species: A) spring-run Chinook Salmon (SRC), B) fall-run Chinook Salmon (FRC), C) Coho Salmon (COS), and D) steelhead (STH). The Augmented Dickey-Fuller test = ADF, Box-Ljung test = ꭓ2, and Breusch-Godfrey test = LM. The horizontal axis shows different lag values, and the vertical axis shows the autocorrelation at each lag value. Also included are 95% confidence limits (bands) for the ACF function correlogram (blue dashed horizontal lines).