Using a mechanistic framework to model the density of an aquatic parasite Ceratonova shasta

The lifecycle of C. shasta

C. shasta is naturally-occurring in river systems across the Pacific Northwest of North America. The parasite alternates between two waterborne spore stages, actinospores and myxospores, that infect salmonid (actinospores) and annelid (myxospores) hosts (Fig. 2; Bartholomew et al., 1997). Actinospores are released into the water from benthic annelids, Manayunkia occidentalis (previously M. speciosa; Atkinson, Bartholomew & Rouse, 2020), encounter the fish host, invade through the gills, move into the circulatory system, and proliferate in intestinal tissue causing enteronecrosis, hemorrhaging, and lesions (Bartholomew, 1998; Bjork & Bartholomew, 2010). Infected fish release myxospores that in turn are consumed by annelids. Actinospores develop within infected annelids and are expelled into the water where they can encounter the next host and complete the lifecycle (Bartholomew et al., 1997). In the Klamath River, actinospores released in the spring can encounter and infect outmigrating juvenile salmon; those released in the fall can encounter and infect adult salmon returning to spawn (Fig. 2). Developmental strategies of C. shasta differ between infections in adult and juvenile salmonid hosts; in adults, release of myxospores occurs after death (i.e., post-spawn), while in smolts and parr the parasite is released shortly after infection (Kent et al., 2014). The parasite released from infected adult carcasses is dispersed downriver where deposition overlaps with habitat for the invertebrate host.

Impacts of C. shasta on salmonids

Fish exposed to high doses of C. shasta are more susceptible to developing clinical disease (Hallett & Bartholomew, 2012; Ray & Bartholomew, 2013; Ray, Holt & Bartholomew, 2012). Juvenile salmonids outmigrating during the spring season in the lower Klamath River can be exposed to high parasite densities. Actinospore density fluctuates temporally and spatially based on the lifecycle described above (see Fig. 2). As river temperatures warm in the spring season, spore densities increase from undetectable levels to peak levels (see Fig. 3 using 2014 as an example; see Bartholomew et al., 2020 for spore monitoring since 2006). The timing of this increase in spores (the “onset” of spores) impacts infection and disease risk for juvenile salmonids and varies as a function of their migration timing (Kovach et al., 2013) through the infectious zone relative to seasonal variation in spore densities (Perry et al., 2019). The highest densities in the year typically occur in the spring (the “peak” spore density; Bartholomew et al., 2020). For example, migration through the infectious zone of juvenile Chinook Salmon of natural-origin spans March–June (David, Gough & Pinnix, 2018); in some years >80% of the run has outmigrated before the onset of spring spores, thereby minimizing disease risk (Som et al., 2019). In other years, outmigration timing coincides with increases of C. shasta densities above levels known to cause high mortality (Ray et al., 2014) likely resulting in the high loss of an entire year class of wild production (Voss, True & Foott, 2018). In addition, approximately five million Chinook Salmon smolts are released from Iron Gate Hatchery (located below Iron Gate Dam, Fig. 1) each year between mid-May and early-June, a period that frequently coincides with peak C. shasta densities (Robinson et al., 2020), exposing hatchery-origin fish to potentially high mortality.

The impacts of C. shasta on outmigrating salmon have been the focus of various management actions on the Klamath River. Federal agencies manage species listed under the Endangered Species Act including Southern Oregon/Northern California Coast (SONCC) Coho Salmon and southern resident killer whales (Orcinus orca) that rely on Chinook Salmon as a primary food source (NMFS, 2019). Flow prescriptions aim to mitigate disease risk to juvenile salmon by disrupting the C. shasta lifecycle (NMFS & USFWS, 2013; NMFS, 2019). Flushing flows (≥171 m³/s for at least 72 h) have been released from Iron Gate Dam to mobilize sediments and disturb river bed stability (Shea, Hetrick & Som, 2016), thereby reducing the abundance of the annelid host (Stocking & Bartholomew, 2007; Alexander et al., 2014; Alexander et al., 2016).

Models developed to simulate components of C. shasta

A model to simulate the timing and abundance of C. shasta in the spring can serve as a critical decision support tool for evaluating the potential effects of flow management scenarios on disease risk of the rearing and outmigration lifestages of juvenile salmon. This model would also be useful in predicting C. shasta-caused mortality in existing management assessments such as the juvenile salmon population dynamics model, the Stream Salmonid Simulator (S3; Perry et al., 2019). Incorporating a simulation model for C. shasta density into S3 would improve model outputs that currently rely on historical observations of spore density that are insensitive to environmental inputs. Simulating spore density patterns that relate directly to disease risk for outmigrating juvenile salmonids is necessary to both evaluate management scenarios and predict their effects on disease-caused mortality of juvenile salmon. The model needs to include (1) the timing of the initial release of actinospores (“onset”) from the annelid host into the water column, (2) the overall seasonal pattern of spore densities, and (3) the peak magnitude of densities for each year (Figs. 3A, 3B and 3C, respectively).

Models have been developed to quantify some phases of the C. shasta lifecycle, providing the foundation upon which to develop a spore density model. A degree-day model using accumulated thermal units (ATUs) predicted the development time of spore stages in hosts (Chiaramonte, 2013; Ray et al., 2015), but this model was based on laboratory experiments and has not been evaluated for accuracy in the development time of C. shasta spores in the Klamath River. The degree-day approach is based upon the temperature-dependency of spore production and release and is commonly used in other myxozoan species to characterize spore development in hosts (e.g., Myxobolus cerebralis, Kerans, Stevens & Lemmon, 2005; De Buron et al., 2017; Tetracapsuloides bryosalmonae, Sudhagar, Kumar & El-Matbouli, 2020). Further, an epidemiological model describing how C. shasta dynamics respond to future climate scenarios found that high winter discharge had the greatest impact on lowering the number of infected juvenile fish (Ray et al., 2015), and a statistical model for annelid habitat suitability related hydraulic and substrate variables to annelid host distribution in the Klamath River (Alexander et al., 2016). Other analyses have found a correlative relationship between the prevalence of infection of hatchery-released smolts and the spring peak of spore densities the following year (Robinson et al., 2020), suggesting seasonal links contributing to spore density. This suite of mechanistic and statistical models provides the groundwork to inform simulations of C. shasta during the outmigration period.

Understanding the bio-physical linkages in the C. shasta lifecycle is important for simulating spore density. A previous attempt to model the abundance of spores in the infectious zone relied on a static population of annelid hosts in a limited section of the river that released actinospores at a constant rate, and the model predicted the prevalence of infection in adult salmon hosts (Schakau, Hilker & Lewis, 2019). However, studies of annelid host populations have documented their range throughout the Klamath River (Stocking & Bartholomew, 2007; Alexander et al., 2014), and laboratory studies have demonstrated spore release rates as highly variable and dependent on temperature and other factors (Bjork & Bartholomew, 2009). In addition, the prevalence of infection in adult salmon is high and not variable, as returning adult fish are immunosuppressed and vulnerable to infection (Dolan et al., 2016); yet a small percentage of carcasses are the highest contributors of myxospores (Foott et al., 2016). As the Schakau, Hilker & Lewis (2019) model predicts infection in adult salmonids, and there is no detectable relationship that links returning spawners with the abundance of spring actinospores (Robinson et al., 2020), a different model design is required to simulate the density of spores in spring that contributes to disease risk in outmigrating juvenile salmonids.

Study objectives

Our objective was to develop a model to simulate the density of waterborne C. shasta during the spring when salmonid outmigration occurs. We aimed to model three key characteristics of spore density patterns that relate directly to disease risk for outmigrating juvenile salmonids. Those characteristics include (1) the timing of the initial release of actinospores (“onset”) from the annelid host into the water column, (2) the seasonal pattern of spore densities, and (3) the peak magnitude of densities for each year (Figs. 3A, 3B and 3C, respectively). Model development was underpinned by data from over a decade of monitoring studies on the Klamath River to explore which biological and environmental components affected these key spore density characteristics. Using these relationships and integrating established mechanistic and statistical relationships, we constructed a model to simulate the density of C. shasta during juvenile salmon outmigration that is responsive to changes in environmental conditions subject to resource management decisions.

Monitoring data compiled for 2006–2018

Multiple monitoring programs collect biological and environmental data on the mainstem of the Klamath River, each with different sampling frequency and methodology described in subsections below (‘Biological variables, ‘Environmental variables’; Table 1). Variables were carefully selected based on the hypothesized potential to influence the lifecycle of C. shasta.

Module A: simulating spore onset using a bio-physical degree-day model

Simulating the initial release date of actinospores each spring relied on coupling parameterized attributes of the parasite lifecycle with monitoring data (salmonid escapement, water temperature, and the occurrence of high flow events) and results from laboratory studies (spore development within annelids). Laboratory experiments demonstrated that following the requisite ATU for spore development (830 ATU; Meaders & Hendrickson, 2009; Alexander et al., 2015), and despite the detection of mature actinospores in infected annelids, release occurred only after temperatures increased above 10 °C (J. Alexander, 2019, unpublished data). To implement these laboratory findings in the module to simulate spore release, both degree-day accumulations and water temperature were included. Given that Klamath River annelids become infected during the fall or winter with myxospores released by the carcasses of adult salmon, in constructing this module, the degree day accumulations were initiated at the beginning of the week when escapement reached the 25th percentile for each year. The 25th percentile of escapement approximates when spawning is underway, the first spawners have died, and myxospore deposition has begun. In the model we assumed after >830 ATU were accumulated and average weekly water temperature ≥11 °C, spores were potentially detectable in water samples; yet for years in which flushing flows occurred, spore release was delayed for 3 weeks after water flow receded below the threshold for flushing flows (Bureau of Reclamation BOR, 2018; Som et al., 2019). Hence, the simulated week of spore release was determined by the timing of adult escapement and spore development, when water temperatures reached at least 11 °C, and whether a flushing flow occurred for 2007–2018. Spawning occurs in the previous calendar year so model output for 2006 was not available. An assessment of the goodness-of-fit of module A compared the simulated week of spore onset to the observed week of spore release from monitoring data. We used a correlation coefficient to assess the strength and direction of the relationship for this and all subsequent goodness-of-fit assessments. We calculated a Pearson correlation coefficient (r) using the cor.test function in R Statistical Software (method = “pearson”; R Core Team, 2018). This comparison excluded 2007 and 2009, when monitoring of spore density started after the onset of spores and therefore the observed timing of release was not available.

Module B: simulating spore density using a temperature-dependent quadratic model

A parabolic thermal response has been described for spore densities based on monitoring studies whereby spore density increases with temperature up to a threshold beyond which spore density declines with further temperature increases (see Perry et al., 2019). The temperature dependency of parasites has been described in other parasite systems with invertebrate hosts by a quadratic relationship between water temperature and abundance (Mangal, Paterson & Fenton, 2008; Paull, LaFonte & Johnson, 2012; Mordecai et al., 2013; Mordecai et al., 2017). Temperature-dependent processes in the lifecycle of C. shasta include the development of spores within hosts (see module A) and also the viability of spore stages in the water once released (Bjork & Bartholomew, 2010). To quantify the independent thermal response isolated from confounding factors, we excluded years in which flushing flow events occurred, as removing hosts by scouring annelids disrupts this relationship, and those in which there was no signal of spore density (selected years: 2008, 2010, 2014, and 2015). The occurrence of flushing flows and low spore density are addressed in module C. The isolated thermal response is important to describe the overall pattern of spore density increase and decrease in response to temperature. Spore densities were normalized to a common relative scale (0 to 1) by dividing spore densities in each year by the annual peak density. A quadratic regression model was then fit to the normalized spore densities with corresponding daily temperature and temperature squared as explanatory variables for the selected years. Using normalized spore densities in the model focuses on describing the shape of the thermal response. Graphical diagnostics were used to evaluate the assumptions of the model; parameters were estimated using the lm function in R (R Core Team, 2018). A goodness-of-fit assessment for module B compared simulated to observed spore densities for years not included in the quadratic regression model. The estimates of the quadratic model were applied to daily temperatures of each year starting at the onset (module A) to simulate normalized spore density.

Module C: simulating peak spore density using a generalized linear model

To simulate peak spore density, we used a generalized linear model (GLM) appropriate for the distribution of spore density data and built upon a published relationship that correlated peak spore density with the prevalence of infection of hatchery-origin Chinook Salmon smolts (HPOI, Robinson et al., 2020). The module framework is described below, followed by model selection and fit.

Spore density data were modeled as generalized linear models with a logarithm link function (Faraway, 2006) to account for the observed mean–variance relationship, whereby variability in spore density increases with mean spore density. We used the negative binomial distribution to account for overdispersion that was evident from fitting Poisson GLMs (Zuur et al., 2009). We assume that peak spore density each year (Y_i; where i represents year) arises from a negative binomial distribution (NB) with mean parameter (µ_i) and dispersion parameter (k): ${\displaystyle Y_i\sim NB\left({\mu }_i,k\right)}$

The effect of covariates on spore density was modeled though a logarithmic function that links a linear function of parameters and explanatory variables to the mean of the distribution: ${\displaystyle log\left({\mu }_i\right)=X\beta }$ where X is a design matrix containing explanatory variable values and β is a vector including an intercept and regression parameter values. We used maximum likelihood techniques to fit models and estimate parameter values using the glm.nb function from the R Statistical Software MASS library (R Core Team, 2018).

Building upon the previously described relationship between peak spore density and the prevalence of infection of hatchery-origin Chinook Salmon smolts (HPOI, Robinson et al., 2020), we evaluated other potential biological and environmental covariates (Table 1) to further explain variation in peak spore density in the spring. The selected variables were limited to those supported by ecological hypotheses. Due to the relatively small number of years available for this analysis, we tested the improvement in model fit provided by adding each selected variable with HPOI against the previous model containing only HPOI. This was done using likelihood ratio tests, which are appropriate for testing among nested models (Lewis, Butler & Gilbert, 2011). The likelihood ratio test is a more conservative choice by providing a formal statistical test of each variable’s improvement in model fit rather than model selection based on the number of units from an Akaike Information Criterion (AIC)-type model selection procedure (see Fig. 1 in Lewis, Butler & Gilbert, 2011). It was decided, a priori, that no models including interaction effects were considered to avoid over parameterizing models given the limited data availability. For each year, peak spore density was simulated by the best model fit. Module C was assessed by a goodness-of-fit comparison between the simulated and observed peak spore density from monitoring data for each year.

Spore density model

The spore density model combines output from modules A, B, and C detailed above. Our goal was to simulate the density of C. shasta spores in the infectious zone through the period of outmigration for juvenile salmonids. The spore density model generates daily values that are averaged for each week 1–30 in 2007–2018; the model terminates at week 30 when most juvenile Chinook Salmon outmigrants have passed downstream of the infectious zone (David et al., 2017), and simulation results for 2006 are not available as the model requires inputs from the previous year. For each year the model combines output from modules A, B, and C. First, spore density is set to zero each week at the beginning of the year before the onset of spore release. The week of spore release is simulated by the bio-physical degree-day module (A). After onset, normalized spore densities (range from 0 to 1) are simulated by the quadratic regression module (B). This module provides the overall seasonal spore density pattern, while the magnitude of spore density for years in which flushing flows occur is simulated by the linear module (C). Last, simulated spore densities are calculated by multiplying the simulated normalized spore density by the simulated peak spore density. A goodness-of-fit comparison between simulated and observed values used the simulated densities for days in which spore density was measured. Simulated spore density was averaged by week.

Module A: simulating spore onset using a bio-physical degree-day model

The simulated onset of spore release in the spring closely matched observations of initial spore detection from monitoring data (Fig. 5A; r = 0.733). For each year, the module started on the average week when escapement reached the 25th percentile, which was week 45 (SD = 0.6). The accumulation of 830 ATU (spore maturity) was reached on average on week 10 of the following year (SD = 2.0). In all cases, this occurred before spores were detected in water samples, which suggests that factor(s) in addition to accumulated temperature influence the timing of spore release. When we included the minimum temperature threshold and delay due to flushing flows, the simulated week of spore release aligned with the observed week of spore release. For 2011, in which observed spore densities remained near zero for the entire spring, the model simulated spore release earlier than the observed week when spores reached 1 spore/L (Fig. 5A), however the model simulated a relatively late week of release. In other years of available data, the week of simulated spore release and observed spore release were similar (mean difference = 0.20 weeks, SD = 1.3).

Module B: simulating spore density using a temperature-dependent quadratic model

The quadratic regression model of spore density and water temperature was based on selected years that did not experience a flushing flow event and that detected spore densities (>1 spore/L; Fig. 6). A quadratic curve fit the non-linear relationship between normalized spore density and water temperature, revealing strong evidence (P < 0.001) that both temperature variables were associated with seasonal spore density (Table 2). There was no evidence of deviation from the assumptions of the model. Goodness-of-fit model assessment compared simulated and observed spore densities (normalized) for years excluded from the development of the quadratic model (Fig. 5B; r = 0.623). Applying the quadratic relationship to all years, daily water temperatures were used to simulate normalized spore density (ranging from 0 to 1) for each day in the spring starting with the week of simulated spore onset (output from module A), and averaged by week.

Module C: simulating peak spore density using a generalized linear model

We analyzed improvements in model fit with the addition of biological and environmental variables that potentially influence the lifecycle of C. shasta (Table 1). Monitoring data for annelid density and prevalence of infection were available for fewer years than spore densities at the time of this study, which would limit the power of analyses. Therefore, these data were not included in subsequent model fit analyses yet should be considered in future years when more monitoring data are available. Potential variables represent the spring conditions of discharge, water temperature, precipitation, reservoir level, and flushing flows. After implementing likelihood-ratio tests and the resulting χ² statistics, a model containing the prevalence of infection of hatchery-origin fish from the previous year (HPOI) with the occurrence of flushing flows (FF) was selected. For most other variables considered, there was no evidence of improved model fit (Table 3). There was weak-to-moderate evidence that annual precipitation improved model fit (P = 0.065), but stronger evidence (P = 0.020) that the flushing flow events improved model fit (Table 3). Given these results, the model containing HPOI and the occurrence of flushing flows successfully simulated the magnitude of peak density for each spring (Fig. 5C; r = 0.776).

The χ² goodness-of-fit test for this model did not indicate any evidence of lack of fit (P = 0.120). For this model, there was strong evidence that HPOI (P < 0.001) and the occurrence of a flushing flow (P = 0.003) were associated with the spring peak of spore densities (Table 4). These results are also consistent with the ecological hypotheses suggesting consideration of these variables in this modeling exercise, as higher HPOI values are correlated with increases in annual peak spore densities, and years with flushing flows are simulated to result in lower annual peak spore densities (Table 4). After accounting for whether years had a flushing flow or not, every 10% increase in HPOI is estimated to increase peak spore density by 102% (95% CI [48–175]% increase). After accounting for HPOI, years with flushing flows are simulated to have peak spore densities that are 90% lower (95% CI [56–98]% lower) than years without flushing flows.

Spore density model

The model of C. shasta generally captured the onset of spore release, seasonal pattern of spore density, and inter-annual variation (Fig. 7), with sufficient contrast to inform management decisions. Most importantly, it closely matches observed values from the onset of spore release and initial increase in spore density for years in which density—and therefore disease risk for juvenile salmonids—was high. The model over- and under-predicted the lowest and highest (respectively) peak spore density values observed (Fig. 5D; r = 0.770), which is common for regression models (Faraway, 2006). For example, the model slightly overestimates spore density in 2011–2013 however both simulated and observed spore densities in these years were <10 spores/L (r = 0.390–0.437) indicating that both simulated and observed spore densities suggest lower disease risk for these years. For years in which spore densities were >10 spores/L, the goodness-of-fit for simulated values compared to observed densities was higher (r = 0.561–0.885).