We use Monte Carlo methods which draw parameters from Bayesian posterior distributions to generate a distribution of population size estimates and trajectories, thus giving managers a fuller accounting of the uncertainty in the population status. We then propagate this population estimate and its associated uncertainty into a model using Monte Carlo methods to assess the impact of fishing bycatch on the species. We show that the population is below the recovery goal of 3,000 adults. The current total population estimate (including juveniles) is approximately 10,000 fish. Our model finds that fishing bycatch pressure reduces an otherwise assumed stable population by a median value of 0.4% per year, which could impede the recovery of the species. Fisheries bycatch is only one of many threats this population faces, and future work is needed to assess how other threats, such as spawning habitat alteration through dams and water diversions, may affect this population’s trajectory. The framework presented here is suitable for further data integration or modular expansion to incorporate the cumulative effects of challenges facing green sturgeon recovery.

We assessed the number of spawning green sturgeon adults in annual surveys of the Sacramento River, CA, USA from the Irvine Finch Boat Ramp (river kilometer 320, just west of Chico) up to Redding (river kilometer 480) (Fig. 1). Acoustic tag data and egg mat studies have confirmed that this is the extent of the spawning grounds (Poytress et al. 2013; Thomas et al. 2014). We surveyed any site along that section of the river with depths greater than 5 m (Erickson et al. 2002). Generally, we see green sturgeon in approximately 40 sites. This section of the river provides the vast majority (effectively all) of the southern DPS green sturgeon spawning locations. 

2.3 Spawner survey

A detailed description of the methods is published by Mora et al. (2015) and is only briefly described here. The survey has taken place continuously since 2010. There are three phases of the survey conducted over three weeks in mid-June. In phase one, a survey crew drifts downstream over the deepest parts of the channel with a depth sounder. The crew contour maps areas of the river with depths greater than 5 m using a sonar system. The survey generally finds approximately 70 areas with a depth over 5 m within the study area (Fig. 1). The crew marks locations with spawner observations in the last 5 years for an automatic revisit in phase three. In phase two, the crew uses a Dual frequency IDentification SONar (DIDSON; Sound Metrics, Bellevue, Washington) video camera to scan for green sturgeon during 3 passes at sites without spawners in the previous 5 years. In phase three, the crew visits all sites where sturgeon have been seen in the past 5 years as well as any new ones added during phase two. At each site, the crew makes 7 passes recording DIDSON footage. The DIDSON footage from phase three is reviewed in random order three times and counts are combined into an estimate of the number of sturgeon at each site location. The sum of counts from all sites is the total number of spawners observed.

2.4 Life table, IPM, and sensitivity analysis

2.4.1 Literature parameters

Both the IPM and life table models need some parameters describing demographics, behavior, and physiology. We took a subset of these parameters directly from the literature (Appendix S1). All these parameters are for the northern DPS green sturgeon as similar data is unavailable for the southern population.

3.4.2 Length vs. age relationship

We used the age vs length data for southern DPS green sturgeon on the Sacramento River from supplement 1 of Ulaski and Quist (2021). We fit these data with Bayesian regression in R using JAGS (packages used rjags, purr, ggplot, dplyr, patchwork, viridis, minpack.lm, ggextra, mcmcplots, and furrr) (Plummer 2003; R Core Team 2015; RStudio Team 2015; Wickham 2016b; Elzhov et al. 2016; Wickham 2016a; Curtis 2018; Wickham et al. 2018; Garnier 2018; Attali & Baker 2019; Pedersen 2019; Vaughan & Dancho 2021). The model and priors are as follows: 

where L is the length, M_L is the mean of the length distribution, T is the precision of the length distribution, L_∞ is the asymptotic length, k is the growth coefficient, A is age in years, and t₀ is agee at zero length. Priors are loosely informed by data from northern green sturgeon (Adair et al. 1983; Farr et al. 2002). We ran three MCMC chains with 1000 adaptation steps and 20000 burn-in steps and saved 10000 samples per chain at 90% thinning. Each chain had random starting values based on draws from the prior distributions. All chains converged based on visual inspection or running means. We compared these results to fits for the northern DPS of green sturgeon as previous population estimates used that data (Beamesderfer et al. 2007; Mora et al. 2018).

3.4.3 Annual survival

Appendix S1 presents an estimate of mortality based on a catch curve analysis from a fishery in the Columbia Estuary. We used the telemetry data from the sturgeon on our system to calculate annual mortality. We then used these two values to bookend the estimate of annual mortality in the life table model and IPM. We used data from the BARD for all green sturgeon tags from 2007-2018. We only used data where length was labeled as either total length or fork length. We converted all lengths to fork lengths and all lengths reported in this manuscript are fork lengths. We used the mean parameters from Eq. 1 to convert the lengths to ages. We grouped the data in 5-year bins to reduce noise. Instantaneous mortality is equal to the slope of the change in counts with age after natural-log transformation. We calculated the slope of the descending arm and converted it from instantaneous mortality to annual survival using (annual survival) = exp(-instantaneous mortality) (Ricker 1975). Subsequent calculations involving survival drew mortality from a uniform distribution over the range between the Columbia River and this estimate.

3.4.4 Probability of being an adult

Rather than using published ages of maturity or the maturity curve implied by the Beamesderfert al. (2007) cohort model, we based the timing of maturity on data specific to the southern population. We calculated the probability that fish of a certain length are adults (i.e. potential spawners) using the same raw data set from the BARD. We flagged fish as adults if they were marked as “mature”, “adult,” or “eggs” (meaning they were caught with eggs) or if detections were above river kilometer 320 (the bottom of the spawning ground) (note the BARD uses a different river kilometer 0, thus river kilometer 320 equates to 410 in the BARD). Only tagging-year records were used so that length and maturity data were contemporaneous. We grouped individuals by sex into females and others (males and unknown). We obtained two separate estimates of the probability of maturation with length using Bayesian logistic regression. The first estimate was a sex-specific hierarchical model fit divided between females and others in which the sexes shared a common slope but had separate intercepts. We used the female parameters from this model in the sensitivity analysis calculation because that analysis needed fecundity. The second estimate was a fit with all the data for use in the population estimate, which needed the total number of spawners. The models are as follows:

where P is the probability of being an adult, M_A is the mean of the probability distribution, L is the length, b₀ is the intercept and b₁ is the slope. Thus, in the hierarchical model, females and others share the b₁ term but have separate b₀ terms. We used broad normal priors (Eq. 2). We ran three chains with 1000 adaptation steps and 20000 burn-in steps and saved 10000 samples at 90% thinning per chain. Each chain had random starting values based on draws from the prior distributions. All chains converged. We compared these results to fits for the northern population of green sturgeon as previous population estimates used that data (Beamesderfer et al. 2007; Mora et al. 2018).

3.4.5 Spawning interval distribution

The population estimate portion of the IPM needs the distribution of spawning intervals for all adults and the sensitivity analysis requires the average spawning interval of females. To calculate these data, we took all detections from the BARD above river kilometer 320 for months during the spawning season (March - September). We removed any detections from the same year the fish was tagged as well as any detections for fish without a detected outmigration between upstream records. For each fish, we then found both the interval between tagging and the first return and, if available, the interval between the first and second return. We divided these into two groups (all fish and females). We then constructed a distribution showing the fraction of fish that have a return interval greater than each interval and calculated the average return interval. We compared these results to fits for the northern population of green sturgeon as previous population estimates used that data (Beamesderfer et al. 2007; Mora et al. 2018).

3.4.6 Calculate the fraction of the population that is adults and spawners

We then sampled from the estimated distributions of survival, maturation, and spawning interval described above to make 10,000 life tables (the average parameter values converged at approximately 8,000 samples) from which we calculated the fraction of the population that is adults. We selected a random set of parameters from the posterior distribution of the Bayesian fit of the length vs. age data, a slope parameter for the length-to-weight equation using a normal distribution and confidence interval from Appendix S1, a mortality value from a uniform distribution between the two annual mortality values, and a random set of parameters from the posterior distribution of the Bayesian fit of the probability of being an adult vs. length data. We used the length age relationship to calculate the length at each age from a life table starting at age 1. Each year we calculated the number surviving using the mortality estimate. We calculated the probability of being an adult given length and thus the total number of adults each year. Using the average spawning interval, we then calculated the number of adults that spawn each year. 

3.4.7 IPM

To calculate 10,000 different rolling population estimates, we first selected values for the number of spawners seen each year. We get these values by drawing a random number for each year of the survey. These random draws are from normal distributions with mean and confidence intervals for the number of spawners observed each year in the survey. We then calculated a rolling estimate of the number of adults using the following equation:

Here N_A is the number of adults in the population, I is the return interval, I_max is the maximum return interval, N_s,i is the number of spawners seen in year j, S is the annual survival rate drawn from a uniform distribution over the estimated range of survival, and P(I > j) is the probability of a return interval greater than j. Years before the first year of the survey (2010) are assumed to have the same values as the first year for calculation purposes. To estimate the distribution of the total population size across years while reflecting parameter uncertainty, we divided each of the estimates of the adult population by the fraction of the population that is adults (mature). We took this fraction from a randomly selected life table drawn from the set of life tables with equal values of annual survival as used in Eq. 3.

3.4.8 sensitivity analysis

To look at the potential effect of fishing bycatch pressure on the population, we first built a set of life tables with a stable population, then applied fishing pressure and calculated the effect on the population growth rate. We took 500 samples from the set of population estimates. Each sample was the population estimate of the most recent year (2018) and the full set of parameters used to calculate that estimate. We then made a life table with the 2018 population estimate and those parameters. We solved for the annual recruitment (egg to first-year survival) to make the life table stable in time. 

To perturb this life table with fishing pressure we used the fleet-wide bycatch estimate for the numbers of southern DPS green sturgeon caught in the open access California halibut fisheries (79% of captures are in this fishery), limited entry California halibut fisheries, and the limited entry and individual fishing quota bottom trawl (Richerson et al. 2019), at-sea hake fishery, and the total population estimate to calculate the possible fractions of the population which are caught. This method allows the number caught to decrease as the population decreases, i.e., we assumed that the bycatch would represent a consistent fraction of the population caught and not simply a consistent number of fish caught regardless of population level. There were only 16 annual bycatch numbers available and so the calculation randomly selected one value each year rather than sampling a fitted distribution. From the bycatch reports we knew the size distribution of bycatch, which we used to distribute the bycatch risk over appropriate length classes. We then simulated each of the 500 sample sets out over 100 years. For each year we drew a random fraction caught from the data set. Then we used that fraction and that year’s current population level to calculate the number caught. We distributed these catches over the lengths in the life table based on the bycatch length data. Each length class got several random draws equal to the total number of fish caught that year. This value was multiplied by the probability of a bycatch incident being in that length class times the probability of the fish dying. The probability of a fish dying is 25% based on the probability of a fish dying in the net (25/508) (Lee et al. 2020 see Figure 2-12 raw data in Appendix S2 of this manuscript) and the probability of a fish dying post-release (11/51) (Doukakis et al. 2020) (1-[1-25/508]*[1-11/51] = 0.25). After the program removed the number of sturgeon bycaught from the population, it advanced all age classes by one year and calculated a new age-1 population. We repeated this process 100 times to project the population over 100 years.