A Spatiotemporal History of Key Pelagic Fish Species in the San Francisco Estuary, CA

Estuaries across the globe have been subject to extensive abiotic and biotic changes and are often monitored to track trends in species abundance. The San Francisco Estuary has been deeply altered by anthropogenic factors, which is reflected in substantial declines in some native and introduced fishes. To track trends in fish abundance, a multitude of monitoring programs have conducted regular fish surveys, some dating back to the late 1950s. While these surveys are all designed to track population-scale changes in fish abundance, they are methodologically distinct, with different target species, varying spatial coverage and sampling frequency, and different gear types. To compensate for individual survey limitations, we modeled pelagic fish distributions with integrated data from many sampling programs. We fit binomial generalized linear mixed models with spatial and spatiotemporal random effects to map annual trends in the spatially explicit detection probabilities of striped bass, Delta smelt, longfin smelt, threadfin shad, and American shad for the years 1980 to 2017. Overall, detection probability has declined by approximately 50% for striped bass and is now near zero for the two smelt species, while threadfin shad and American shad have both experienced fluctuations with only slightly reduced detection probabilities by 2017. Detection probabilities decreased dramatically for these fishes in the Central and South Delta, especially after the year 2000. In contrast, Suisun Marsh and the North Delta acted as refuge habitats with reduced levels of decline or even increased detection probabilities for some species. Our modeling approach, using disparate datasets, demonstrates the simultaneous spatially driven decline of pelagic fish species in a highly altered estuary.


Introduction
Estuaries are highly productive and often urbanized systems located at the interface between freshwater and marine environments. Biological productivity is fueled through the input of both terrestrial and marine nutrients, as well as through increased residence time due to salinity-driven density gradients and tidal forcing (Pérez-Ruzafa et al. 2011). Estuaries support diverse assemblages of aquatic species along the salinity gradient, from obligate stenohaline species at the marine and freshwater fringes to oligohaline species that can utilize the entirety of the estuary (Whitfield et al. 2022). The high productivity and location of estuaries at the interface between freshwater and marine environments has also resulted in the general colonization of estuaries by humans seeking the benefit of rich food sources and protected ports (Wilson 1988;Lotze et al. 2006;Cabral et al. 2022).
The dynamic nature of estuaries has encouraged, and in some cases necessitated, environmental modification for human use. Many large estuaries have been diked, drained, and dredged for urbanization, transport, water management, flood control, and agriculture. Inflows have been diverted or impounded in reservoirs (Cabral et al. 2022). These changes have substantially altered abiotic conditions and have often changed the habitat in estuaries to favor non-native species (Cabral et al. 2022; Moyle and Stompe 2022.). Species Communicated by Mark S. Peterson 1 3 introductions are common in estuaries through ballast water exchange from international shipping traffic, as well as by the deliberate introduction of recreational, ornamental, and bait/forage species (Moyle and Stompe 2022).
Estuaries are frequently monitored to track trends in the abundance and distribution of aquatic species of recreational, commercial, and cultural importance, many of which have declined because of extensive abiotic and biotic changes (Blaber et al. 2022;Cowley et al. 2022). Estuarine monitoring is undertaken by diverse state, federal, tribal, academic, and/or non-governmental organizations, all with potentially different objectives (Anderson 2005). The oftenpiecemeal implementation of surveys complicates traditional analyses of species trends. In some cases, this results in underutilization of data due to concerns about differences in methodology and bias (Stompe et al. 2020;Huntsman et al. 2022). Unfortunately, logistical and/or financial limitations of surveys means that trends in abundance and distribution of estuarine species are often incomplete when based on analysis of a single survey.
The San Francisco Estuary (Estuary) is a system with many independent, long-term fish monitoring programs. The Estuary includes the tidally influenced portions of the Sacramento and San Joaquin Rivers, the Delta, Suisun Bay and Marsh, San Pablo Bay, Central and Southern San Francisco Bay; it terminates at the Golden Gate (Fig. 1). The Estuary has been highly altered by myriad abiotic and biotic changes since large-scale colonization by European-Americans in the 1800s. Abiotic changes include extensive physical and hydrologic alterations through water diversions, levees, floodplain and marshland reclamation, hydraulic mining, and dams. Dams are now present on every major river in the Estuary's 163,000 km 2 watershed; they have severely disrupted natural flow and sediment regimes (Cloern and Jassby 2012; Whipple et al. 2012;Schoellhamer et al. 2013;Herbold et al. 2014). The Estuary has also been subject to numerous species invasions through ballast water dumping in its freshwater and saltwater ports, and statesponsored and illicit intentional introductions (Cohen andCarlton 1998, Moyle andStompe 2022). Introduced species include fishes that have altered predation and competition dynamics (Grossman 2016;Moyle et al. 2016;Nobriga and Smith 2020), invertebrates that have converted productivity pathways (Greene et al. 2011), and plants that have altered turbidity and water velocity throughout much of the Estuary . As a result of these changes, and other stressors such as overfishing (Yoshiyama et al. 1998) and pollution (Brooks et al. 2012), many native and some introduced fish species have experienced drastic declines in the past 100 years, and especially the last 40 years Stompe et al. 2020).
Along with numerous smaller diversions, the Estuary contains two major water export facilities located in the South Delta, the State Water Project (SWP), and Central Fig. 1 Simplified spatial plane of the San Francisco Estuary with applied barrier components. White background is wetted area and grey background is land. Blue dots represent the center of "water" spatial mesh triangles and green dots represent the center of "land" spatial mesh triangles. Select cities surrounding the Estuary, the location of the Montezuma Slough salinity control gates, the location of the Delta Cross Channel, and the location of the South Delta export facilities (State Water Project, Central Valley Project) are included for reference. Numbered regions are identified for regional descriptions of distribution trends. 1 = South San Francisco Bay, 2 = Central San Francisco Bay, 3 = San Pablo Bay, 4 = Carquinez Strait, 5 = Suisun Marsh (top) and Suisun Bay (bottom), 6 = Sacramento-San Joaquin Confluence, 7 = North Delta, 8 = Central Delta, 9 = South Delta 1 3 Valley Project (CVP; Fig. 1), which export a large proportion of freshwater inflow for agricultural and municipal use (Gartrell et al. 2017;Moyle et al. 2018). To track the impacts of water infrastructure operations on Estuary fish species, state and federal agencies and a research group at the University of California, Davis (UC Davis), have operated regular fish surveys dating back to the late 1950s (Honey et al. 2004;Herrgesell 2012;Tempel et al. 2021). Early surveys primarily focused on tracking the abundance of young of year striped bass, an introduced yet recreationally and culturally important species within the Estuary (Turner and Chadwick 1972;Chadwick et al. 1977;Stevens et al. 1985). Later, additional surveys were started to track the abundance and survival of outmigrating juvenile salmonids (Dekar et al. 2013) and to track changes in fish and invertebrate assemblages (Baxter et al. 1999;Matern et al. 2002;O'Rear et al. 2021) rather than single species.
While most Estuary fish surveys are designed to track fish abundance over a wide spatial expanse, they are logistically and economically restricted in total number of stations and frequency of sampling events. In addition, surveys sample different micro-and macro-habitats due to differences in gear type, spatial coverage, and project goals (Tempel et al. 2021). As the micro-and macro-distribution of fishes shift in response to abiotic and biotic change, so does the ability of any given survey to capture the fishes (Sommer et al. 2011). Due to these differences in survey coverage and fish distribution, multiple surveys and gear types must be analyzed in tandem to fully understand trends in Estuary fish abundance and distribution.
In this paper, we leverage an integrated dataset of eight long-term surveys (Stompe et al. 2020) to examine trends in the distribution and probability of detection of five important pelagic fish species: striped bass (Morone saxatilis), Delta smelt (Hypomesus transpacificus), longfin smelt (Spirinchus thaleichthys), threadfin shad (Dorosoma petenense), and American shad (Alosa sapidissima). We use spatially explicit species distribution modeling to fit binomial generalized linear mixed models to the integrated survey data and then (1) identify key areas of importance for these species over time, (2) pinpoint time periods of major shifts in distribution and abundance, and (3) describe the effects of freshwater outflow on distribution.

Study Species
Five species were chosen for our analysis (striped bass, Delta smelt, longfin smelt, threadfin shad, American shad) due to their ecological, cultural, and recreational significance in the Estuary. These species represent important recreational fisheries (striped bass and American shad) and both native (Delta smelt, longfin smelt) and introduced (threadfin shad) forage fishes. In addition, all these fishes require productive estuarine pelagic environments during part or all their lives (Moyle 2002), so their abundances can be indicators of the "health" of pelagic habitats.
Striped bass is an introduced, relatively long-lived, semianadromous, and fecund species that relies on productive estuaries for rearing (Raney 1952). Since their introduction (1879; Dill and Cardone 1997), they have been one of the primary catch species in agency and university surveys, although catches have declined considerably over the years (Kohlhorst 1999;Sommer et al. 2007). Many of the surveys included in our analysis were established specifically for the capture of juvenile striped bass.
Delta smelt is a small osmerid endemic to the Estuary. They are generally an annual species and are an obligate estuarine species (Moyle et al. 1992;Moyle 2002). Despite once high abundance, they are now rarely caught by Estuary surveys and are listed as threatened under the Federal Endangered Species Act (FESA) and endangered under the California Endangered Species Act (CESA; Tempel et al. 2021).
Like Delta smelt, longfin smelt are small native osmerids; they are found along the Pacific Coast of North America, are more halophilic than Delta smelt, and can live 2 to 3 years (Moyle 2002). Historically, they were highly abundant within the Estuary but have since declined and are now relatively rare . As a result, they were listed as threatened under the CESA in 2009 (Tempel et al. 2021).
Threadfin shad are introduced, small, deep bodied clupeids that typically live two to three years (Moyle 2002). Despite their somewhat recent introduction to the system (1962;Feyrer et al. 2009), threadfin shad have also experienced declines in abundance .
American shad are another introduced clupeid species (introduced 1871; Dill and Cardone 1997), but they reach larger sizes (adults 30-60 cm FL) than threadfin shad and generally migrate into the Pacific Ocean after rearing in the Estuary (Moyle 2002). Previous work has not identified major reductions in Estuary American shad abundance and instead has seen recent increases in angler catch (Ferguson 2016).

Survey Data
The surveys included in our modeling effort are the Califor-  Table 1). Of the stations we include in our analysis, there is considerable spatial overlap between the surveys in the San Pablo, Carquinez, Suisun, and Sacramento-San Joaquin Confluence Regions ( Fig. 1 -regions 3-6; Fig. S1). Conversely, the Bay Study Otter Trawl and Bay Study Midwater Trawl are the only surveys with stations in the Central and South San Francisco Bays ( Fig. 1 -regions 1 and 2; Fig. S1) and only the Fall Midwater Trawl, Summer Townet Survey, Beach Seine Survey, and Chipps Island Trawl have stations in the Delta ( Fig. 1 -regions 7-9; Fig. S1). The longest running of these surveys (Summer Townet Survey) started in 1959 and the most recent (Bay Study Otter and Midwater Trawls) in 1980, and all have operated continuously on at least an annual basis through 2017. Several of these surveys were originally designed to describe and track trends in young of year striped bass abundance (Fall Midwater Trawl, Summer Townet Survey), some were designed to track the outmigration of juvenile salmonids through the Estuary (Chipps Island Trawl, Beach Seine Survey), and others were designed to track assemblages of fish and invertebrate populations (Suisun Marsh Otter Trawl/Beach Seine, Bay Study Otter, and Midwater Trawls). All were implemented to determine the effects of water diversion and/or entrainment into export facilities on fish populations. These surveys primarily capture small and/or juvenile fishes due to their specific gear types and netting mesh sizes. For this reason, surveys are generally able to describe trends in the relative abundance of small fishes such as Delta smelt, longfin smelt, and threadfin shad, but only represent juvenile trends of large fishes, such as striped bass and American shad.
We integrated data from these surveys into an aggregate dataset for the years 1980 through 2017, retaining key variables such as date, coordinates, and number of individual fish captured by species (Stompe et al. 2020). For consistency of annual spatial extent, we only include those survey stations which were sampled at least once annually in our analyses. In addition, several Beach Seine Survey stations located on the Sacramento River upstream of the Delta were omitted because of their distance from other downstream stations. Because effort shifted among certain years, data used from the Chipps Island Trawl was seasonally restricted to April through June and from the Fall Midwater Trawl to September through December to standardize seasonal effort. All other surveys generally operated year-round.
Catch per unit effort was calculated as total number of fish caught per seine or trawl, with a total of 103,341 sampling events. While other Estuary data integration efforts have instead chosen to index catch by the volume of water sampled (Huntsman et al. 2022), not all surveys we include record this metric, nor is it as meaningful of a metric for some gear types (i.e., beach seines). In addition to indexing catch per trawl or seine, we included species presence or absence for each sample. The inclusion of a binary metric allows for the modeling of probability of detection, regardless of water volume sampled.
In this integrated format, these data represent approximately 40 years of trends in Estuary fish abundance at a much greater seasonal and spatial density than could be provided by any single survey. In addition, the breadth of gear types included in the eight-survey dataset mean that benthic, pelagic, and littoral species are all targeted to some degree.

Data Analysis
Using the aggregate eight-survey dataset, we modeled spatiotemporal trends in the probability of detection of the five pelagic fish species using the package "sdmTMB" (Anderson et al. 2022). We constructed and fit binomial generalized linear mixed models (GLMMs) of species presence by maximum likelihood for each fish species from the aggregate dataset across a restricted spatial mesh. Model predictions from the binomial models represent the probability of detection rather than the total predicted catch of a given species. Because of this, these models do not necessarily capture absolute changes in abundance. However, binomial model predictions do provide information on changes in abundance because survey gear is more likely to detect species at higher densities assuming that they are relatively evenly distributed within the sampled habitats and assuming limited density dependence of catch, such as gear saturation at high fish densities (Godø et al. 1999). During model construction, Tweedie distributions (Tweedie 1984) were also tested, but model fit was poor with this distributional family (Hartig 2022).
A spatial mesh for efficiently modeling spatial and spatiotemporal autocorrelation was generated using the "cutoff" method, with a minimum of 2-km spacing between mesh vertices ("knots") and resulting in a mesh with 179 knots (Anderson et al. 2019). Knot spacing was iteratively chosen to reduce model overfitting in highly sampled areas of the Estuary (Suisun Bay, Carquinez, etc.), while also providing acceptable spatial coverage in less intensively sampled areas. The spatial mesh was geographically restricted to the wetted area of the Estuary by restricting spatial autocorrelation between geographically close, yet ecologically distinct, habitats ( Fig. 1). Shorelines were simplified using barrier polygons and small channels were expanded to allow the mesh to fit an adequate number of knots for even representation of sparsely sampled areas. Likewise, only major islands separating heterogeneous habitats were included to increase the number of knots in the sinuous regions of the Delta and Suisun Marsh. Finally, the Montezuma Slough Salinity Control Gates and Delta Cross Channel Gates were treated as open to reflect the average condition during sampling periods (Fig. 1). Simplification of the barrier polygon was an iterative process, as early barrier polygons with full shoreline and channel complexity had very few knots present in the Delta and Suisun Marsh regions.
The notational structure (Eq. 1) of the binomial GLMMs is shown below. The number of sampling events that detected a given species at location s in year t and month m by sampling program p, y s,t,m,p , is modeled as a binomial random variable with expected probability of detection at a single sampling event equal to s,t,m,p and N s,t,m,p sampling events over the month. We used a logit link to model s,t,m,p as a function of independent year effects ( t ), survey effects ( p ), and a cubic spline for month (s(m)). The variable s represents spatial random effects, s,t represents spatiotemporal random effects, and t represents the spatially varying coefficients through time Y (scaled year t centered around zero with a standard deviation equal to one). The random effects ( s , s,t ) and the spatially varying coefficient ( t ) are drawn from Gaussian Markov random fields with Matérn covariance matrices Σ , Σ , and Σ , respectively (Barnett et al. 2021). Spatial random effects were included to account for unmeasured variables that are approximately fixed through time (depth, distance upstream, substrate, etc.) whereas spatiotemporal random effects were included to account for unmeasured variables that are likely to change over both time and space (salinity, temperature, food availability, etc.; Anderson et al. 2022). Spatiotemporal random effects were treated as independent and identically distributed.
Residuals from the non-random effects were then simulated by drawing 500 samples from the fitted model and tested using the "DHARMa" package in R (Hartig 2022). Residual (1) y s,t,m,p ∼ Binomial N s,t,m,p , s,t,m,p , logit s,t,m,p uniformity was tested using the one-sample Kolmogorov-Smirnov test, residual dispersion was tested using the DHARMa nonparametric dispersion test via the standard deviation of residuals fitted versus simulated, and residual outliers were tested using the DHARMa outlier test based on the exact binomial test with approximate expectations (Hartig 2022). None of the models demonstrated unacceptable levels of residual uniformity, dispersion, or outliers, indicating good fit.
Previous publications have identified the potential pitfalls of generating models using disparate datasets in an integrated format (Walker et al. 2017;Moriarty et al. 2020;Huntsman et al. 2022). When unaccounted for, differences in survey effort, gear efficiency, and overall catchability can introduce significant biases in abundance and spatiotemporal density trends (Walker et al. 2017;Huntsman et al. 2022). However, our inclusion of survey as a fixed effect accounts for these biases, allowing separate intercepts to be fit for each of the eight surveys. In addition, we test for potentially different catch trends between pelagic and littoral or shallow sampling gear types by plotting smoothed trends (generalized additive model; Wickham 2016) in CPUE (as catch per trawl/seine pull) of several geographically close stations from the Beach Seine Survey, Suisun Marsh Fish Study, Fall Midwater Trawl, and Summer Townet Survey. Species chosen for comparison of CPUE trends were selected based on regional results from model outputs.
Once models were fit, we made predictions of the probability of detection and spatially explicit slopes of predictions over time for each of the five species across a 500-m grid of the wetted area for the visualization of Estuary wide trends. We then generated estimates of the spatial slope standard deviation as a metric of uncertainty by sampling (n = 200) from the joint precision matrix.
Once prediction dataframes were made, we calculated mean estimates of detection probability by decade for each spatial point and rasterized means into smooth prediction planes using the R package "ggplot2" (Wickham 2016). We plotted rasterized spatial slopes, representing relative change through time, along with the standard deviation of the estimates of the spatial slopes. We assumed a linear relationship between detection probability and time for the spatial slopes, a potential simplification of trends which does not identify step-changes but can be compared to decadal estimates to confirm validity. We also plotted mean annual estimates of the probability of detection at grid points for each of the five species by applying a smoothing function (generalized additive model) by year.
To measure distributional sensitivity to changes in Delta outflow conditions (amount of water exiting the Delta after water exports) and overall shifts in population distributions, we calculated the annual predicted center of gravity (COG) and 95% confidence intervals along a longitudinal gradient 1 3 for each of the modeled fish species. Delta outflow data was sourced from the California Department of Water Resources (CDWR 2022) as millions of acre-feet of water per calendar year. COG is a metric which represents the mean location of a population, weighted by density of observations, and although it is imperfect at describing local trends and/or detecting changes at distributional extremes (Barnett et al. 2021), it can be a useful metric for measuring population movement along a distributional gradient (Thorson et al. 2016). Due to computational limitations, COG estimates were generated using predictions at survey station points rather than at all points on the prediction grid. As a result, the COG of each species does not necessarily represent the true longitudinal center of each population, but changes in COG over time and relative differences between species are valid given the temporal consistency in spatial sampling. We calculated COG longitudinally to best reflect the flow direction of the Estuary, which generally runs East-West from the Delta through the Central San Francisco Bay (Fig. 1).
To test for effects of Delta outflow on COG, temporal trends in COG, and differences in COG between species, we constructed a generalized additive model using the package "mgcv" (Pseudo-R Code, Eq. 2; Wood 2011). The point estimates of COG for each species from the sdmTMB prediction outputs were included as the response variable, with yearly estimates weighted by their variance. We applied smooth functions with thin plate basis splines (bs = "tp") to "Year" by "Species" and "Delta Outflow" by "Species," and "Species" was included as a linear fixed effect. Finally, we printed generalized additive model results in tabular format and generated plots of yearly COG point estimates, 95% confidence intervals around the point estimates, the fit trendline and 95% confidence interval by species from the generalized additive model, and the spline effect of Delta outflow in million-acre feet on COG (Wickham 2016;Coretta 2022).

Model Results
All models fully converged and fit the data acceptably well as determined through residual testing in the DHARMa package (Hartig 2022). The results of residual dispersion tests and residual outlier tests were significant for some models (Table S1); however, this was driven by the exceptionally large number of data points included in the model rather than poor model fit (Hartig 2022). For example, the dispersion value of 0.992 for the American shad model was significant (p = 0.012) as was the outlier test (p = 0.008) with just 467 outliers out of 103,341 observations. The results of model diagnostics are included in the supplementary material (Fig. S2, Table S1).
The output of the GLMMs showed differing coefficients by "Survey" for all species as well as negative coefficients for the linear component of the smooth effect of "Month" sampled for threadfin shad and American shad (Table 2).

3
The differing coefficients by "Survey" indicate differences in the catchability of species by survey methodology and gear type, while the negative coefficient of the linear component of "Month" is likely due to seasonal differences in either local abundance, absolute abundance, or vulnerability of capture across life stages for the two shad species. As expected, the pelagic gear types used by surveys such as the Bay Study Midwater Trawl, the Fall Midwater Trawl, and the Chipps Island Trawl were generally most effective at catching the pelagic species included in our analysis as indicated by model coefficients, with some notable exceptions (Table 2). For example, the Fall Midwater Trawl, a survey specifically designed to capture young of year striped bass, had a lower model coefficient than the reference survey (Bay Study Midwater Trawl). Beach seines and benthic trawls typically had negative coefficients among the modeled fish species, indicating lower capture efficiencies than the reference survey (Bay Study Midwater Trawl). Despite this, the Beach Seine Survey and Suisun Marsh Otter Trawl both showed similar trends in catch as was seen for geographically close pelagic trawls (Fig. S3).

Spatial Trends
In general, species show an overall reduction in spatial probabilities of detection over the modeled time period (Fig. 2). Spatial slopes are negative in most regions for most species and are exclusively negative or near zero for Delta smelt and longfin smelt (Fig. 3). Positive slope values are present in limited regions for striped bass (Suisun Marsh, North Delta, South San Francisco Bay), threadfin shad (Confluence, Suisun Bay, Suisun Marsh), and American shad (North Delta, Suisun Bay, Suisun Marsh, South San Francisco Bay).
Species show some shared patterns in distributional changes in the probability of detection over time; most notable of which is the reduction in estimates in the Central and South Delta ( Fig. 1 -regions 8-9; Fig. 2). This trend exists for striped bass, threadfin shad, and American shad which historically had relatively high probabilities of detection in these regions, but which had very low detection probabilities by the 2010s. Further supporting this are the spatial slopes which are strongly negative in the Central and South Delta for these species (Fig. 3). The reduction in detection probability in these regions drives a constriction in distribution away from large parts of the Delta and towards the Confluence and Suisun regions ( Fig. 1 -region 6 and 5).
Unlike the other three species, the two smelt species were rarely detected in the Central and South Delta regions at any point during the modeled time period. Longfin smelt were mostly distributed downstream, with historically high probabilities of detection in the Central San Francisco Bay through the Confluence region ( Fig. 1-regions 2-5). Conversely, Delta smelt were not found in the lower portions of the Estuary and were instead most likely to be detected in the North Delta, Confluence, and Suisun Regions (Fig. 1 -regions 5-7). Both species seem to exhibit a reduction in overall detection probabilities rather than a constriction of distribution.
Another notable trend is the persistently higher probability of detection in the Suisun Marsh and Suisun Bay regions ( Fig. 1 -region 5 top and bottom, respectively) relative to other parts of the Estuary for all species. For most species, it appears that the probability of detection does not increase in Suisun over the decades, but rather it decreases less. This is supported by the spatial slope plots which show relatively less negative slope in the Suisun Bay and especially Suisun Marsh regions (Fig. 3). American shad was the one species which did have positive slopes throughout much of the Suisun Marsh region, indicating an increased detection probability between 1980 and 2017.
In general, there was little uncertainty in the spatial slopes in highly sampled regions where species were detected by the eight surveys over the modeled time period. Uncertainty was high at the edge of the spatial mesh, such as the north end of San Pablo Bay ( Fig. 1 -region 3), or where species were never or rarely detected (Delta smelt, Central and South San Francisco Bays, Figs. 2 and 3). Spatial slopes should be interpreted cautiously in these specific areas given the relatively high level of uncertainty.

Detection Trends
The overall probability of detection for the included species generally declined between 1980 and 2017 (Fig. 4), evidence that abundance may be declining. Declines are most evident for striped bass and the two smelt species, the latter of which are now rarely detected in the eight surveys. Conversely, the detection probabilities for threadfin shad and American shad appear to somewhat rebound near 2017 after lows in 2010 through 2012.
The trends in detection probability are primarily nonlinear, with intermittent periods of increase or stabilization. Striped bass have the largest overall decline, from a detection probability of approximately 0.35 in the early 1980s to less than 0.15 by 2012 (Fig. 4). Delta smelt have an initial steep decline in the early 1980s, followed by relatively stable to increasing detection probability until another period of steep decline in the early 2000s. After this decline, Delta smelt again stabilize until the mid-2010s, at which point they decline to near zero. Longfin smelt trends are similar to Delta smelt, but with a less dramatic initial decline followed by a low in the early 1990s and a rebound in the late 1990s before ultimately declining to near zero as well. Threadfin shad trends are somewhat similar to longfin smelt, but as stated earlier, they have partially recovered in the years since 2010. Finally, American shad are unique in their trends, with a somewhat stable detection probability until a decline between 2005 and 2010, followed by a recovery.

Center of Gravity
The generalized additive model (Eq. 2) of the smooth effects of year and Delta outflow by species on COG, and of the linear effects of species on COG, fit with an adjusted R-squared of 0.969 and explained 97.7% of the deviance in the data. The smooth term for Delta outflow indicated strong non-linear effects of Delta outflow on COG for Delta smelt (edf = 5.10, F = 7.77, p < 0.001), longfin smelt (edf = 6.92, F = 6.09, p < 0.001), and threadfin shad (edf = 7.68, F = 9.29, p < 0.001). Conversely, effects of Delta outflow on COG were linear and weak or non-existent for striped bass (edf = 1.00, F = 3.39, p = 0.068) and American shad (edf = 1.00, F = 2.42, p = 0.122).
The COG of each fish species partitioned by easting, a Universal Transverse Mercator (UTM) coordinate system measurement roughly equivalent to longitude. Threadfin shad were distributed the furthest upstream (Table 3, Fig. 5), followed by Delta smelt, American shad, and striped bass clustered within approximately 5 km of one another, and longfin smelt the furthest downstream (Table 3, Fig. 5).
The point estimates and modeled fit of COG by species also resulted in several distinct patterns over the modeled time period. The smooth term for year indicated strong non-linear effects of year on COG (Table 4) for longfin smelt (edf = 8.52, F = 13.66, p < 0.001) and threadfin shad (edf = 8.61, F = 13.20, p < 0.001), and weak, semi-linear effects for Delta smelt (edf = 2.88, F = 2.42, p = 0.069). There was little effect of year on COG and a linear or nearly linear trend for striped bass (edf = 1.00, F = 0.51, p = 0.475) and American shad (edf = 1.78, F = 1.99, p = 0.131). Longfin smelt had the most dramatic temporal trends in COG, spanning approximately 20 km over the modeled time period (Fig. 5). In addition, longfin smelt COG remained further downstream during the period after 2002, despite some periods of extreme drought (CDEC 2021). Threadfin shad also showed strong temporal trends in COG, but with a unique pattern where the population generally shifted upstream during the period from approximately 1990 to 2010 (Fig. 5). Finally, Delta smelt remained relatively centered within a 5-km band, with little interannual variability (Fig. 5).

Discussion
Using new developments in spatiotemporal modeling, we leveraged the rich but fragmented monitoring data in the Estuary to demonstrate changes in the spatial and temporal  probability of detecting key pelagic fish species during a period of considerable abiotic and biotic change. Large swaths of the Estuary that historically supported high detection probabilities of striped bass, threadfin shad, and American shad, including the South and Central Delta, are now relatively devoid of these species, driving population constrictions to the Suisun and Confluence regions. Over the same time period, Delta smelt and longfin smelt experienced relatively even declines throughout the Estuary.
The detection probability of all species declined to some extent from their levels in 1980; however, the trends and future outlooks differ by species. The differences in life history and age-structured distribution of species color the interpretation of our results. Striped bass are a long-lived semi-anadromous species most likely to be caught by survey gear during their first year of life; results reflect juvenile distribution and are not directly indicative of adult behavior or abundance. Likewise, American shad typically migrate out of the Estuary (and into the Pacific Ocean) after their first year (Carothers et al. 2021), so juveniles are best represented in our results. Longfin smelt also sometimes leave the Estuary; however, they are susceptible to survey gear throughout their lives so results may be interpreted as representing the total local population. Next, threadfin shad remain vulnerable to survey gear throughout their short (1-3) year lifespans so survey catch is representative of population trends, although their population may also be supplemented by individuals that are flushed into the Estuary from upstream reservoirs during wet years. Finally, Delta smelt are a small annual species (Moyle 2002) that are vulnerable to survey gear throughout their lives and are fully restricted to the sampled estuarine areas, so our results may be interpreted as representing the spatiotemporal trends of the species as a whole.
Given these species-specific differences in model interpretation, it is clear that Delta smelt and longfin smelt have declined precipitously since the 1980s (Fig. 4). These species are now rarely detected by Estuary surveys. The threadfin shad population has also experienced declines in overall detection probability; however, it appears to be more robust in its ability to shift to different regions as evidenced by positive spatial slopes in the Confluence and Suisun regions (Fig. 3). As a result, threadfin shad overall probability of  detection has somewhat recovered since 2010 (Fig. 4). The most dramatic decline in detection probability is seen for striped bass (Fig. 4), indicating a reduction in either spawning success or juvenile survival over the modeled time period. American shad spawning success or juvenile survival has also somewhat declined between 1980 and 2017, although, their probability of detection is now only slightly below historic highs after a recovery since 2010 (Fig. 4).
The lower level of overall decline in juvenile American shad versus juvenile striped bass, despite similar distributional patterns, may be somewhat driven by increased utilization of Suisun Marsh by American shad as the Central and South Delta became inhospitable (Fig. 3). Sensitivity of annual COG to outflow conditions indicate different relative effects of climate and water management for each of the species. Delta outflow has major effects on the location of the salinity gradient, and thus the highly productive low-salinity zone (MacWilliams et al. 2015). Given that longfin smelt COG are associated with Delta Outflow (Table 4, Fig. 5), this indicates that the distribution of this species may annually shift to track areas of favorable salinity and/ or productivity. Conversely, species such as American shad whose COG are insensitive to different outflow conditions (Table 4, Fig. 5) may not be as plastic in their annual distribution, potentially due to reliance on fixed habitat features rather than water conditions. It is difficult to identify from our analyses whether longitudinal plasticity in response to changes in Delta outflow is advantageous, as species with highly variable COG, such as longfin smelt, and species with relatively stable COG, such as striped bass, have both experienced dramatic declines in their probability of detection. However, previous work has identified outflow as positively associated with abundance for longfin smelt and to a lesser degree American shad and striped bass (Colombano et al. 2022).
Trends in annual COG over time suggest differential responses to changing environmental conditions and highlight life history differences among the species (Table 3, Fig. 5). The highly variable annual COG for longfin smelt and threadfin shad support a plastic response in distribution, or potentially differences in success between multiple subpopulations within the Estuary. For example, multiple spawning populations of longfin smelt have been identified within the Estuary (Lewis et al. 2020), so regional differences in spawning success or survival could shift annual COG (Colombano et al. 2022). Conversely, the relatively fixed annual COG for Delta smelt, striped bass, and American shad indicate a general reliance on fixed habitat features as well as either a lack of subpopulation structure or subpopulations with similar interannual spawning success or survival.
The abiotic and biotic drivers of the described changes in detection probability and distribution are likely complex and interacting. For example, over the modeled time period, the Estuary saw changes in water export regimes (Gartrell et al. 2017), the introduction of several highly invasive plant and invertebrate species (Cohen and Carlton 1998), recordsetting droughts , and extremely wet years (CDEC 2021). These factors interact, changing the amount and quality of habitat for native and introduced pelagic fishes.
For example, invasive plants such as Brazilian waterweed (Egeria densa) have benefitted from reduced turbidity due to upstream impoundments and the constant freshwater condition maintained by water export operations . Dense stands of Brazilian waterweed have reduced water velocity in some areas, dropping out additional suspended particulate matter and capturing nutrients from upstream sources (Yarrow et al. 2009;Durand et al. 2016). This has resulted in potentially reduced pelagic productivity (Vanderstukken et al. 2011;Durand et al. 2016), a shift in zooplankton communities important for small and larval fish diets (Espinosa-Rodríguez et al. 2021), and reduced turbidity (Hestir et al. 2016). These factors make small fishes more susceptible to predation (Ferrari et al. 2014) and have been shown to negatively affect both striped bass and Delta smelt ). There are myriad examples of such interacting and cascading effects that have reduced the suitability of the pelagic habitat within the Estuary (Brown and Moyle 2005;Sommer et al. 2007;Brooks et al. 2012;Cloern and Jassby 2012;Sabal et al. 2016).
An overarching trend in the distribution and regional abundance of these species is the relative insulation of the Suisun Region, and to a lesser degree, the North Delta, from overall declines in detection probability (Figs. 2, 3). There are many potential drivers of this, one of which is likely the historically lower levels of colonization of submersed aquatic vegetation, such as Brazilian waterweed, in these regions. Brazilian waterweed is largely limited by salinity in Suisun Marsh and Bay (Borgnis and Boyer 2016) and was previously limited by high turbidity in the North Delta . Given the relatively lower levels of fish decline in detection probability in these regions, they could prove important for maintaining viable populations of pelagic fishes in the future.
Another potential driver of the insulation of the Suisun and North Delta regions may be due to changes in habitat use and the relative capture efficiency of the different surveys. Sommer et al. (2011) identified a shift of age zero striped bass away from deep channel habitats to shallower shoal habitats through analysis of the Fall Midwater Trawl and Bay Study Midwater Trawl data. This shift indicated a potential behavioral response, possibly due to decreased pelagic food availability, but did not fully explain the marked decline in striped bass catch. Given our inclusion of surveys and gear types that mainly or exclusively sample shallow habitats, a behavioral response may indeed explain some of the trends in detection probability that we identify.
In our analysis, Suisun Marsh was primarily sampled by the UCD Suisun Marsh Fish Study and the North Delta by the USFWS Beach Seine Survey (Fig. S1). On average, these surveys sample more shallow habitats than do the other surveys due to the use of smaller vessels (Suisun Marsh Fish Study) and depth limited gear (Beach Seine Survey). This could potentially explain some of the reduced levels of decline in these regions, provided that the species have shifted their distribution from channels to shoals, small sloughs, or other shallow habitats. While we acknowledge this as one potential driver of our results, similar trends in catch of geographically close Beach Seine Survey, Summer Townet Survey, Suisun Marsh Fish Study, and Fall Midwater Trawl stations (Fig. S3) indicate that sampling methodology is not solely responsible for the trends observed in our model output. This is further supported by Mahardja et al. (2017), which also found decreasing CPUE of American shad and Delta Smelt (although not meeting the author's threshold for significance for Delta smelt) among Beach Seine Survey stations in the Delta.
Despite lower levels of decline or even increased detection probability in the Suisun region and the North Delta, these regions may not necessarily represent ideal habitat in the face of system wide degradation. These regions may simply be better than those regions which have experienced substantial declines in detection probability (Central Delta, South Delta). Under this scenario, fish may be shunted away from previously productive habitats into regions which have experienced relatively less change. This likely partially explains the increased detection probability of some species in Suisun Bay and Marsh and the North Delta over the modeled time period.
We also acknowledge that regional trends in detection probability were sometimes non-linear (Fig. 2), complicating the interpretation of spatial slopes. For example, American shad increased in detection probability in the 1990s and 2000s before declining in the 2010s, resulting in a positive spatial slope in this region (Fig. 3). Given this, the spatial slopes are a useful yet imperfect tool for describing regional changes in detection probabilities. Spatial slopes show overall changes in regional detection probabilities between 1980 and 2017, but do not identify periods of major step-change. For this reason, mean decadal detection probabilities (Fig. 2) are more useful for identifying time periods of change and should be interpreted along with spatial slopes to fully understand the trends in detection probability.
Station density is somewhat sparse in the North and South Delta, with one or few surveys representing catch in these regions. Specifically, trends in the North Delta are most influenced by catch from the USFWS Beach Seine Survey, while trends in the south Delta are most influenced by catch from the Summer Townet Survey (Fig. S1). These regional differences in station density and representation should be considered when interpreting the absolute detection probability of species such as Delta smelt in the North Delta, although annual CPUE at paired Beach Seine Survey and Summer Townet Sites ( Figure S3) indicates similar trends between surveys.

Conclusions
By leveraging existing long-term survey data in an integrated modeling format, we have described trends in the distribution and detection probability of five pelagic fish species in the Estuary. The modeling techniques we employ have most commonly been used to describe trends in large adult marine fishes, but we demonstrate their ability to model trends in juvenile or small estuarine fishes as well. Our approach also demonstrates the value of using an integrated data set due to the greatly increased spatial density and coverage. These data can detect distributional trends that would otherwise not be covered by a single survey. We are aware that measures must be taken to ensure that modeling with disparate data does not impart unacceptable biases, so we included "survey" as a fixed effect and only used consistently surveyed stations. These measures should be sufficient to control for the methodological disparities. Our modeling and data integration methods should not only prove useful for management of Estuary fishes, but also for describing trends in distribution and abundance of fishes in other inland, estuarine, and marine systems with multiple independent surveys.
The individual long-term fish surveys of the Estuary have collected valuable data for tracking trends in the distribution and abundance of the species we considered here, and indeed for most fish species found within the Estuary. While any one survey can describe part of a species' story, it is only when surveys are analyzed in concert that we can see the true extent of change. The increased spatial breadth and detail of an integrated analysis allows us to see much more granular and localized changes in distribution.
The Estuary has experienced dramatic changes to its hydrology, biotic communities, and physical structure, which in turn has reduced the detection probability and distributional breadth of both native and naturalized pelagic fish species. The species analyzed here include fishes that hold considerable ecological, recreational, and cultural value among California stakeholders. We show, through distributional shifts and spatial slopes, that a major driver in the reduced overall detection probability of pelagic fishes is the decline in their probability of detection in large portions of the Delta. Conversely, Suisun Marsh and the North Delta appear to function as refuge habitats for at least a few of the species, reinforcing that they should be managed as high priority refuges for fish conservation.
While our analyses identify regions and time periods of change for these fish species, they do not specifically 1 3 identify biological or abiotic drivers of these changes. In future efforts, our models may be expanded through the inclusion of spatially explicit data for both biotic and abiotic predictors, including bathymetry, temperature, and LIDAR imagery of submersed aquatic vegetation. Expansion of our models in this way would further refine our understanding of important habitat criteria for pelagic fishes in the Estuary, which should result in more effective prioritizing of habitat restoration and conservation efforts.
Acknowledgements Caroline Newell and Dr. Eric Ward provided helpful reviews of an early version of this manuscript, as did two anonymous reviewers. Avery Kruger contributed to the original idea for our analysis. This study would not have been possible without the numerous staff and managers of Estuary fish surveys who spent considerable time and effort collecting the data that we used.
Funding This project was funded in part by the California Department of Fish and Wildlife and the California Department of Water Resources.

Data Availability
The data used in this study are available by request from the corresponding author (DKS).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.