Long-term trends in seasonality and abundance of three key 1 zooplankters in the upper San Francisco Estuary 2

(300 words)


Introduction
Zooplankton are critical components of aquatic food webs, connecting primary producers to upper trophic species such as fishes.Most fishes rely on zooplankton food in their larval stages when the starvation risk is highest (Hunter 1981).The seasonality of fish reproduction is often timed such that larvae coincide with peaks in the abundance of their zooplankton prey (Cushing 1969;Cushing 1990).However, when the seasonality of zooplankton abundance shifts due to climate change, species introductions, or other factors, larval starvation risk can increase as their peak abundance no longer coincides with peak food availability (Edwards and Richardson 2004;Durant et al. 2007).Thus, an understanding of the patterns in zooplankton abundance is critical to understanding the drivers of fish abundance.
The San Francisco Estuary (estuary) is home to several fish species listed under the United States Endangered Species Act and/or the California Endangered Species Act.Most depend on zooplankton prey for part or all of their life cycle.One of these species, Delta Smelt (Hypomesus transpacificus) is endemic to the estuary and relies on zooplankton for its full life cycle (Brown et al. 2016).Another fish, Longfin Smelt (Spirinchus thaleichthys), also relies on zooplankton throughout its life cycle (P.Chigbu and Sibley 1998;Paulinus Chigbu and Sibley 1998;Barros et al. 2022).A number of fish species, including Delta Smelt and Longfin Smelt, declined sharply in abundance in the early 2000s, during the "pelagic organism decline" (Feyrer et al. 2007;Sommer et al. 2007).This decline is thought to have been caused in part by reduced zooplankton food supply (Winder and Jassby 2011;Brown et al. 2016;Moyle et al. 2016).
Zooplankton research in the estuary to date has focused on a few key taxa, most notably the copepods Pseudodiaptomus forbesi, Eurytemora affinis, and Limnoithona tetraspina, as they are important food for threatened Delta Smelt and Longfin Smelt.(Bouley and Kimmerer 2006;Kimmerer et al. 2014;Kayfetz and Kimmerer 2017;Kimmerer et al. 2018).In floodplains adjacent to the Delta, some cladocerans have also received attention, particularly Daphnia spp., due to their importance in the diet of juvenile salmonids (Goertler et al. 2018;Corline et al. 2021).
One prior study has investigated changes in zooplankton phenology in the estuary.Merz et al. (2016) used a high-level approach to evaluate changes (from 1972-2014) to the date of maximum abundance of five zooplankton taxa: E. affinis, Pseudodiaptomus spp., other calanoids, cyclopoids, and non-copepods.They found that the peak abundances of all these zooplankton groups except the cyclopoids have shifted weeks to months earlier.
However, analyses of the date of peak abundance can overlook the presence of multiple seasonal abundance peaks.Furthermore, the coarse taxonomic resolution of this study excluded analyses of some key species.Many members of the zooplankton community remain understudied, despite their role in the pelagic food web In this study, we extensively review prior knowledge and derive new insights from monitoring data on three key but understudied zooplankton taxa.We chose to focus on Bosmina longirostris, Acartiella sinensis, and Acanthocyclops spp.(Fig. 1) due to their importance in estuarine ecology and fish diets as well as the paucity of prior studies on their dynamics in this estuary.We first review past studies from this estuary and other regions where these taxa occur, then describe the importance of each taxa in fish diets by tabulating their records in fish diets, and lastly examine their long-term trends in seasonality and abundance in the estuary by applying Bayesian generalized additive mixed models to an integrated database of zooplankton monitoring data (Bashevkin et al. 2020;Bashevkin et al. in review).Prior research on the study species Bosmina longirostris B. longirostris (Fig. 1A) is the most abundant cladoceran in the freshwater reaches of the Sacramento San Joaquin Delta (Bashevkin et al. 2020;Jeffres et al. 2020), where it has been a major component of the zooplankton community over the past 40 years (Ambler et al. 1985).Like Daphnia spp., it is abundant in off-channel habitat in the estuary (Corline et al. 2021), but unlike Daphnia spp. it is also abundant in the larger channels of the South Delta (Bashevkin et al. 2020;Jeffres et al. 2020).Bosmina spp. is consumed by fish in the estuary and elsewhere, including juvenile salmonids and Delta Smelt (Roegner et al. 2015;Goertler et al. 2018;Slater et al. 2019), but it rarely occurs in diets of larval Longfin Smelt (Jungbluth et al. 2021;Burris et al. in review), likely due to a mismatch in the salinity preferences of Bosmina spp.and Longfin Smelt.Its small size makes it generally less preferred than larger Daphnia spp. in diets of juvenile salmon (Craddock et al. 1976;Holm and Møller 1984), and diets with high percentages of Bosmina spp. in juvenile perch result in low growth rates (Romare 2000).Both Daphnia spp.and Bosmina spp.tend to have lower nutritional quality (as measured by fatty acids) than copepods (Kratina and Winder 2015).
B. longirostris is a small, filter-feeding cladoceran found in lakes and rivers throughout the world, where it has received little attention due to its small body size (Adamczuk 2016).Despite its small size in comparison to members of the better-studied genus Daphnia, B. longirostris can still have a large impact on the food web.B. longirostris feeds broadly on phytoplankton, protists, and bacteria, ranging from 1-15 µm, and preferentially consumes algae over bacteria, when available (Balcer et al. 1984;Onandia et al. 2015).B. longirostris is particularly effective in depressing biomass of small phytoplankton (Carpenter and Kitchell 1984) and is an efficient consumer of ciliates and bacteria (Jürgens et al. 1996), playing a key role in the microbial loop.It can also thrive on many types of cyanobacteria (Tõnno et al. 2016).B. longirostris is more tolerant of disturbance than many species of Daphnia (Hart 2004;Adamczuk 2016), with greater resistance to toxic cyanobacteria (Matveev and Balseiro 1990;Jiang et al. 2013;Jiang et al. 2014;Jiang et al. 2017), and a higher salinity tolerance than many other cladocerans (Adamczuk 2016).Bosmina spp.also has a higher thermal tolerance than many other freshwater zooplankton (Drenner et al. 1981;Jiang et al. 2014), leading to increasing dominance over other zooplankton as temperatures rise.B. longirostris typically lives 20-50 days and reproduces parthenogenically, producing 2-6 broods of 2-4 embryos each, though temperature, salinity, and predation pressure may impact life span and reproduction (Adamczuk and Mieczan 2019).

Acanthocyclops spp.
Acanthocyclops spp.(Fig. 1B) is a cyclopoid of unknown origin which occurs mostly in freshwater (Orsi and Mecum 1986).Previous studies in the estuary identified the species of Acanthocyclops in this region as Acanthocyclops vernalis; however, research in other areas discovered that A. vernalis is a species complex consisting of three cryptic species; Acanthocyclops robustus, Acanthocyclops americanus, and A. vernalis (Alekseev et al. 2002;Dodson et al. 2003;Miracle et al. 2013;Alekseev 2021).Due to possible morphological misidentifications, it is not known whether all three species in the A. vernalis complex are native to the estuary and were present in the past; however, Jungbluth et al. (2021) did find molecular evidence of all these species in larval Longfin Smelt diets recently collected in this region.
Before the introduction of the small cyclopoid L. tetraspina in 1993 (Orsi and Ohtsuka 1999), Acanthocyclops spp. was the most abundant cyclopoid in the estuary (Orsi and Mecum 1986).Acanthocyclops spp. is an important component of fish diets in this region, particularly Longfin Smelt (CDFW unpublished data; Hobbs et al. 2006;Burris et al. in review) and larval Pacific Herring (Clupea pallasii) (Jungbluth et al. 2021), with Acanthocyclops spp.being the most common cyclopoid consumed by larval Longfin Smelt (Jungbluth et al. 2021).Acanthocyclops spp.also has a higher nutritional value than L. tetraspina, and a similar fatty acid composition to that of the calanoid copepods E. affinis and P. forbesi (Kratina and Winder 2015).
While research on the ecology and biology of Acanthocyclops spp. is relatively limited in the estuary, the A. vernalis complex has been studied extensively in other regions including Europe, Russia, and the Great Lakes.Studies in these areas show that, despite the morphological similarities, A. vernalis, A. robustus, and A. americanus have different ecologies, life cycles, and environmental tolerances (Alekseev et al. 2002;Miracle et al. 2013;Alekseev 2021).A. vernalis and A. robustus inhabit littoral or near-benthic areas and A. vernalis has a benthic naupliar stage (Alekseev et al. 2002;Miracle et al. 2013) and vertically migrates at night (Evans and Stewart 1977).Both A. vernalis and A. robustus are also predatory, consuming copepod nauplii, cladocerans, rotifers, and occasionally larval fish (Anderson 1970;Kerfoot 1978;Gliwicz and Stibor 1993;Piasecki 2000).A. americanus, however, is pelagic throughout its lifecycle, found in varying salinities from the Mediterranean Sea (Alekseev 2021) to freshwater lakes (Alekseev et al. 2002), and is omnivorous.A. americanus nauplii consume primarily algae, with later life stages also consuming filamentous algae and cyanobacteria, in addition to cladocerans, nauplii, and rotifers (Enríquez-García et al. 2013;Sarma et al. 2019).All species can produce about 100 eggs per female, develop to sexual maturity in 10-14 days, and live 30-75 days depending on conditions, with A. americanus growing and reaching maturity faster than the others (Alekseev 2021).Species in the A. vernalis complex likely have different environmental tolerances, as has been shown by studies of their seasonal and spatial variation in abundance in other regions.A. vernalis may be more tolerant of colder temperatures than the other species, whereas A. americanus could have a higher temperature tolerance based on laboratory experiments and timing of peak abundance in areas outside the estuary (Alekseev 2021).

Acartiella sinensis
In the fall of 1993, the non-native calanoid copepod A. sinensis (Fig. 1C,D) was first detected in the estuary (Orsi and Ohtsuka 1999).Likely introduced via the ballast water of ships, the large (~1.2 -1.5 mm in length) predatory calanoid is native to the estuaries and coasts of Southeast Asia (Shen and Lee 1963;Srinui and Ohtsuka 2015).The species has been recorded from estuaries along the East China Sea in salinities around 18-21 ppt (Shen and Lee 1963) to the brackish marshes of Thailand in salinities around 5 ppt (Srinui and Ohtsuka 2015).In the Pearl River estuary of China, A. sinensis was the dominant copepod in brackish waters with salinities less than 15 ppt (Tan et al. 2004).Sampling in the Thale-Noi lake of Thailand showed changes in temperature and salinity were the main environmental variables impacting densities of A. sinensis in the region (Inpang 2008).
Within a year after invasion, A. sinensis became the second most common calanoid copepod in the upper estuary, with its highest abundances in the lower salinity zone of Suisun and the West Delta during summer and fall (Hennessy 2018).The predatory A. sinensis has been shown to feed on the nauplii and copepodid life stages of other copepods in the estuary, primarily those of the abundant L. tetraspina and P. forbesi (York et al. 2014;Slaughter et al. 2016;Kayfetz and Kimmerer 2017).A. sinensis has also become a food item for the endangered Delta Smelt, appearing in their diets in summer and fall (Slater and Baxter 2014;Slater et al. 2019).
The invasion of A. sinensis was not without consequence.Since the invasion in 1993, the zooplankton assemblage in the low-salinity Suisun area has shifted in trophic composition.Once dominated by herbivorous cladocerans and copepods such as E. affinis and P. forbesi, the community has become more "top-heavy" with the spread of A. sinensis (Kratina et al. 2014;Kratina and Winder 2015).The sustainment of A. sinensis in the lowsalinity zone and its high predation rate on the nauplii of P. forbesi is linked to a shift in the distribution of P. forbesi out of the low-salinity zone and upriver into more freshwater habitats (Kayfetz and Kimmerer 2017).This shift in the distribution of the important P. forbesi could have implications for the majority of planktivorous fishes that feed on the calanoid copepod populations.

Diet data
To investigate the occurrence of our study species in fish diets, we summarized data from the California Department of Fish and Wildlife (CDFW) Diet and Condition Study ("Diet Study", https://wildlife.ca.gov/Conservation/Delta/Special-Studies).The CDFW Diet Study identifies and enumerates gut contents of young pelagic fishes in the estuary including Delta Smelt, Longfin Smelt, Striped Bass (Morone saxatilis), Threadfin Shad (Dorosoma petenense), American Shad (Alosa sapidissima), and other species of interest collected by various Interagency Ecological Program monitoring surveys (Slater and Baxter 2014;Hammock et al. 2017;Slater et al. 2019).Fish were weighed, measured, the guts dissected, and prey items were identified to the lowest possible taxonomic level.In some cases, prey items were heavily digested and therefore could only be identified to higher taxonomic levels than our study taxa.Thus, these prey items were not included in our data.Fishes were chosen for processing based on other study needs, so relative sample size for each fish is not correlated to abundance in the estuary, but provides some information on relative rates of zooplankton consumption for each species.
The number of fish processed, number with food in their stomachs, fork length range, and number of fish that had eaten our focal taxa were extracted from the database.To maintain the same spatial coverage as the zooplankton datasets, only fish caught in San Pablo Bay and upstream were included in this analysis.

Zooplankton data
The data used for these analyses were obtained from an integrated database of five long-term zooplankton monitoring surveys in the upper estuary.These include the Environmental Monitoring Program (EMP), 20-mm Survey (20mm), Fall Midwater Trawl (FMWT), Summer Townet (STN), and Fish Restoration Program monitoring (FRP).These surveys are described in detail in Kayfetz et al. (2020) and Bashevkin et al. (in review).
Briefly, EMP samples monthly year-round since 1972, 20mm samples every other week March through July since 1995, STN samples every other week June -August since 2005, FWMT samples monthly September through December since 2011, and FRP samples annually to monthly near tidal marshes (or areas soon to be restored to tidal marshes) March through December since 2015.Each survey samples at a set of fixed stations (Fig. 2).EMP also samples at a set of moving stations (entrapment zone, EZ) at locations where the bottom conductivity is 2 and 6 mS/cm.Many of these surveys collect other parameters such as fish abundance or water quality, but only the time period of zooplankton sampling is described above.Furthermore, sampling locations and frequency have changed over time (Kayfetz et al. 2020;Bashevkin et al. in review).The data from each survey are integrated by the R package zooper (Bashevkin 2021), which imports and standardizes the data from each survey.For the purposes of these analyses, we selected data from the mesozooplankton size-class, which corresponds to samples from the modified Clarke-Bumpus nets used by each survey, with mesh sizes of 150 -160 µm.We selected all available data from each of our study species without missing values in our covariates.Adults and juveniles were analyzed separately since these distinct life stages have different behaviors and abundance drivers and play distinct demographic roles.A. sinensis was introduced to the estuary in 1993 (Orsi and Ohtsuka 1999) Exploratory data visualization revealed tighter relationships of each species' abundance with log-transformed salinity than raw salinity values.Thus, salinity was natural log-transformed for analyses.Furthermore, we centered all covariates (including logtransformed salinity) by their mean and standardized them by their standard deviation.Lastly, since many of the sampling stations from the different surveys were nearby one another (Fig. 2), we clustered all stations into groups within a 1 km radius.

Model structure
We fit Bayesian generalized additive mixed models with the R package brms (Bürkner 2017;Bürkner 2018), which uses the Bayesian modeling platform Stan (Stan Development Team 2021), as well as the R package mgcv (Wood 2011) to handle the smooth construction.Models were fit with a hurdle lognormal distribution, using the catch per unit effort (CPUE; organisms m -3 ) of the specified taxa and life stage as the response variable.A hurdle model was used to account for the large number of 0s in the data.
Our overall approach was to model the probability of presence with a smoothed function of salinity, and the non-zero abundances with smooth functions of day of year, salinity, year, and their interactions, while accounting for space with a random intercept for each station.Our combination of a Bayesian method, which propagates uncertainty and handles unbalanced data, with a generalized additive model approach, which accounts for key covariates like salinity and seasonality, allowed us to resolve inconsistencies in the sampling designs, and incorporate both fixed stations and stations that move with the salinity field.
The hurdle probability (probability of 0 CPUE) was modeled with a cubic regression spline of salinity with a low basis dimension (k) of 5 since the relationship was expected to have a simple shape.The abundance of Acartiella sinensis was so strongly seasonal that we modified the hurdle component to also include seasonality.Thus, for A. sinensis we modeled the hurdle probability with a two-dimensional tensor product smooth of salinity (cubic regression spline, k=5) and day of year (cyclic cubic regression spline, k=4).
The non-zero CPUEs were modeled with a three-dimensional tensor product smooth of day of year (cyclic cubic regression spline, k=13), salinity (cubic regression spline, k=5), and year (cubic regression spline, k=5).The basis dimension for day of year was set to 13 to match the monthly nature of these sampling programs (since the smooth is cyclical, a basis dimension of 13 has 12 independent functions).The basis dimension for salinity was set to the low value of 5 because the relationship with salinity was expected to be a simple unimodal shape, and a basis dimension of 5 would still allow much greater complexity than that.The basis dimension for year was set to the low value of 5 because we were interested in evaluating broad long-term patterns, rather than fine-scaled year-to-year abundance trends.
Thus, the results of this model represent broad long-term trends, not predicted abundances on specific years.For juvenile A. sinensis, the basis dimension for year was reduced to 3 since they have only been counted in samples since 2006 and thus the timeseries is shorter.We also included a random intercept for each station cluster.
We fit separate models on each species and life stage (4 total).Models were run on three chains, each for 5,000 iterations including 1,250 used for the warmup that were then discarded.We used weakly informative priors as recommended by the Stan authors (Stan Development Team 2021).
Each model was validated and checked prior to use.All models were inspected to ensure adequate sampling by verifying the posterior effective sample size (> 100 per chain) and Rhat values (< 1.05) (McElreath 2015).We further inspected visual plots comparing the model outputs to the raw data to ensure they matched.These plots included the proportion of zero values, the distribution of non-zero CPUE values, and the predicted non-zero CPUE values for each row in the dataset.We also inspected the spatiotemporal variograms for spatiotemporal autocorrelation using the R package gstat (Pebesma 2004;Gräler et al. 2016).
We detected some residual autocorrelation, and thus used a conservative significance threshold of p < 0.01 to account for any potential impacts.

Model predictions
We visualized predicted values from the models to explore the abundance patterns of each species and life stage.Predicted values were generated over a range of covariates that included all combinations of 6 evenly spaced time-points per month (from 1972 to 2019) and a series of salinity values selected as quantiles from the raw data (every 0.05 from 0.05 to 0.95).Since there were some gaps in the timeseries (e.g., winters were not sampled some years in the 1970s and 1980s), those same gaps are preserved in the model predictions to avoid extrapolation.Model predictions were then plotted along with their 99% credible intervals.To visualize the multidimensional model outputs, we created three plots for each set of model predictions.Each plot had one of the covariates (salinity, day of year, or year) on the x-axis while the other two variables were illustrated with color or separate plots.For the two covariates included as color or separate plots, we chose a subset of the unique values used in generating model predictions in order to reduce plot complexity and aid interpretation.These values were chosen as an evenly spaced subset of the values available.
To explore spatial patterns in abundance, we extracted the mean estimated value from each station cluster random intercept.These values were then plotted over a map of the study region.

Diet results
The CDFW Diet Study has processed 11,301 fish caught in San Pablo Bay and upstream through the Delta consisting of nine species; Delta Smelt, Longfin Smelt, age-0 Striped Bass, Threadfin Shad, American Shad, Mississippi Silversides (Menidia beryllina), Pacific Herring, Prickly Sculpin (Cottus asper), and Tridentiger spp.gobies (Table 1).A wide range of fork lengths was processed for all species, except for Pacific Herring, Prickly Sculpin, and Tridentiger spp.Gobies which only included larval fish.Time and seasonality of collection varied by fish species with year-round diet information only available for Delta Smelt and Longfin Smelt.
Our study zooplankton were eaten by at least one individual of each fish species, except Pacific Herring, which consumed none and is thus not discussed further here (Table 1).Acanthocyclops spp.were present in the most fish stomachs (n= 885) and A. sinensis juveniles were eaten the least (n= 88).Of the fish that consumed Acanthocyclops spp., Delta Smelt had the highest consumption (n= 800, 21% of fish with food present), and Striped Bass the least (n= 3 fish, 0.3%).Neither Prickly Sculpin nor Tridentiger spp.consumed Acanthocyclops spp.
Bosmina was present in the guts of all fish species (n= 235) and was the third most consumed of our study zooplankton.Threadfin Shad ate Bosmina the most (n= 53, 9%), Delta Smelt the second most (n= 163, 4%), and Longfin Smelt ate the least (n= 1, 0.04%).
Table 1: Diet results for fish processed by the CDFW Diet and Condition Study and the number of fish that consumed each zooplankton species."# fish" represents the number of fish of each species that were processed for stomach contents."Season" indicates the seasons in which fish were collected, with "Sum" = summer and "Win" = winter.Fork lengths are given in mm.Zooplankton species are represented as follows: A = Acanthocyclops spp., AS = A. sinensis adults, AS juv = A. sinensis juveniles, and B = B. longirostris.
Percentages represent the percent of fish with food present in their stomachs that consumed each prey species.

Bosmina longirostris adults
The abundance of adult B. longirostris has regularly peaked in the Spring between April and May.In earlier years, they had another large peak in the late summer to fall between August and October.In some years (1980 and 1985), this fall peak was as large as the spring peak.However, by 1990 the fall peak was greatly reduced to just a small increase, which has since continued to decrease in size, becoming non-existent by 2015 (Fig. 3).
B. longirostris adults are most abundant in freshwater, and abundance decreases as salinity increases (Fig. S1).In the moderate salinity of 1.113, the fall peak was larger than the spring peak in earlier years.However, as with the lower salinities, the fall peak was greatly reduced by 1995, shrinking smaller than the spring peak even as the overall abundance in this salinity continued to decrease over time (Fig. 3, S1).longirostris adult abundance, after controlling for the other covariates (day of year, salinity, and year).

Acanthocyclops spp. adults
The abundance of Acanthocyclops spp.adults has peaked regularly in the spring from April to May in all years.This peak is apparent in the two lower salinity levels (0.062 and 0.137) but generally not in the two higher salinities (1.402 and 13.017).In the second highest salinity of 1.402, abundance peaked in the winter from February through March in most years, although limited sampling in these months in some years may have masked the signal (Fig. 6).
Acanthocyclops spp.were most abundant in the lower salinities, peaking around 0.3.
Abundance decreased on either side of that peak, but fell much lower in the highest salinities (Fig. S2).Over time, the relationship of Acanthocyclops spp. with salinity leveled out such that they became equally abundant in all salinity levels in most months (Fig. S2).
The abundance of Acanthocyclops spp.has dramatically declined over time in all months.There was a slight uptick in the 1980s to 1990s in most months, but populations crashed again after this period (Fig. 7).Overall trends in January and February appear mostly flat with generally low abundance in all years but the missing data in those months may be obscuring patterns.
The highest Acanthocyclops spp.abundance after controlling for the other covariates was in Cache Slough, Suisun Marsh, the Southeastern Delta, Carquinez Strait, San Pablo Bay, and the Napa River.The lower Sacramento and San Joaquin rivers through Suisun Bay generally had lower abundances, as did most of the western-most areas (Fig. 8).for a full description.

Acartiella sinensis adults
A. sinensis had the strongest seasonality of the three species investigated.Adults peak in the fall from August through December with the highest levels in September and October.
Abundances then dip close to zero from February through May (Fig. 9).
Adult A. sinensis were most abundant in moderate salinities between about 1 and 4.
The effect of salinity on abundance has increased over time, especially in May through July where the peak was greatly reduced in earlier years (Fig. S3).
Unlike the other two species, A. sinensis adults did not exhibit any overall long-term decreases in abundance.However, the timeseries is shorter, only starting in 1994.In most months, the most recent abundance is similar to the earliest abundance level, but abundance did increase over time in March through July.The abundance peaked in the 2010s in most months and some months also had an earlier peak around 2000 (Fig. 10).
Spatially, A. sinensis adults had the highest abundance (controlling for all other covariates) along the corridor from the lower Sacramento River just below Cache Slough all the way through to Carquinez Strait.The Southeastern and Northern Delta had the lowest abundances (Fig. 11).Figure 11.Estimated values of the clustered station random intercepts for Acartiella sinensis adults.See Fig. 5 for a full description.

Acartiella sinensis juveniles
Like the adults, juvenile A. sinensis also had strong seasonal abundance patterns, peaking over just a few months and then subsiding to close to zero abundance.Peaks occurred in the summer from July through September, but the width of the seasonal peak grew over time.In 2006 they were abundant for just 2 months (July and August) while in 2015 and 2018 they were abundant from April through November (Fig. 12).
A. sinensis juveniles were abundant in higher salinities > 4 but declined at the very highest salinities close to 16.Their abundance in lower salinities increased over time but always remained lower than the higher salinities (Fig. S4).In most years, the seasonal abundance peak was 1-2 months later at the highest salinity of ~16 than the other salinity levels (Fig. 12).
While the timeseries was much shorter (2006 -2020) for A. sinensis juveniles than any of the other species and life stages investigated, we did detect some long-term trends in a few months.The trends were most apparent in the second highest salinity of 2.703 where they were most abundant.Abundance increased over time in April through June and decreased over time in August.This corresponds to the widening of the seasonal peak over time.The other months generally did not have any significant long-term trends (Fig. 13).
The spatial pattern of A. sinensis juveniles was less clear than the other species and life stages.However, they were generally most abundant along the San Joaquin River corridor in the Southern Delta and in some Suisun Bay stations.They were least abundant in the lower Sacramento River between Cache Slough and the Confluence, as well as in the Napa River and Eastern Suisun Marsh (Fig. 14).
Figure 12.Seasonal patterns in Acartiella sinensis juvenile abundance with 99% credible intervals.See Fig. 3 for a full description.

Discussion
We found marked changes in the seasonality and overall abundance of three key zooplankton species in the upper estuary.B. longirostris no longer peaks in abundance in the fall months, Acanthocyclops spp.dramatically declined in all months and lost its strong relationship with salinity, and A. sinensis adult abundance has become increasingly driven by salinity while juveniles have developed wider abundance peaks.In this process, we have documented the relationship of each species with salinity and seasonality back to the beginning of monitoring or their introduction, increasing our understanding of their ecology and importance in the estuary.

Importance in fish diets
Of our study species, Acanthocyclops spp. was consumed by the most fish, particularly Delta Smelt and Longfin Smelt.This is consistent with other studies that show the importance of Acanthocyclops spp. in Longfin Smelt (Hobbs et al. 2006;Jungbluth et al. 2021), Delta Smelt (Lott 1998;Hammock et al. 2017), and larval Pacific Herring (Jungbluth et al. 2021) diets.Acanthocyclops spp.has a similar fatty acid composition to the calanoid copepod E. affinis (Kratina and Winder 2015), the preferred prey of Delta Smelt (Slater andBaxter 2014), andLongfin Smelt (Barros et al. 2022;Burris et al. in review), and presents an additional prey resource during the critical life stages of these fishes.
Adult A. sinensis were consumed the second most of our study species, but juveniles were consumed the least.Since its introduction, A. sinensis has become important prey for many fishes in the estuary during fall, especially Delta Smelt (Slater and Baxter 2014;Hammock et al. 2017;Slater et al. 2019).Because A. sinensis is predatory, it has a lower nutritional value (i.e. higher Carbon to Nitrogen ratio) than herbivorous copepods such as P.
forbesi (Kratina and Winder 2015), another common prey of Delta Smelt. A. sinensis predates and limits the availability of P. forbesi in the low salinity zone (Kayfetz and Kimmerer 2017), replacing P. forbesi as a prey source for juvenile and young-of-year Delta Smelt. A.
sinensis is not common in Longfin Smelt or Striped Bass diets, likely due to the peak seasonal abundance occurring when both fish species are larger and have transitioned to larger prey such as mysids and shrimp (Feyrer et al. 2003;Barros et al. 2022, CDFW unpublished data).
B. longirostris were consumed the second least of our study species.Threadfin Shad was the most frequent consumer of B. longirostris compared to the other fish species, although the amount was still low despite their spatial overlap.Unlike most of the other fishes examined, Threadfin Shad switch their feeding strategy from filter feeding to benthic organisms and detritus based on prey availability (Ingram and Ziebell 1983) and this could explain the low number of individuals that consumed B. longirostris.B. longirostris, and other cladocerans, have a lower nutritional value than copepods (Kratina and Winder 2015).
However, cladocerans are still an important prey resource for fishes in the freshwater regions in the estuary such as salmonids (Goertler et al. 2018;Aha et al. 2021), Mississippi Silversides, and Threadfin Shad (Whitley and Bollens 2014), or during times when copepods are less abundant such as for Delta Smelt in winter and spring (Hammock et al. 2017;Slater et al. 2019).
Based on previous research (Lott 1998;Hobbs et al. 2006;Hammock et al. 2017;Jungbluth et al. 2021), it is likely that the presence of our study species in fish diets is higher than our data indicates.This could be due to damage during digestion hindering identification of prey to the taxonomic level of our study taxa, or that most fish samples were limited to specific seasons and may not have been collected at a time that our zooplankton species were present in the environment.Thus, our results likely represent undercounts of the actual abundance of our study taxa in fish diets.
Seasonal, salinity-related, and geographic abundance patterns Currently, B. longirostris and Acanthocyclops spp.adults peak in the spring while A.
sinensis adults peak in the fall and juveniles peak in the summer.The spring peaks line up with the spawning of Delta Smelt while summer and fall peaks provide critical food for Delta Smelt juveniles and young-of-the-year (Slater and Baxter 2014;Slater et al. 2019).The spring peaks also correspond to the outmigration of juvenile salmon and could provide important food necessary to increase growth rates and reduce oceanic predation risks (Herbold et al. 2018;Phillis et al. 2018;Zeug et al. 2019).
regions. A. sinensis juveniles had some similarities to the adults but were much more abundant along the San Joaquin River corridor and had low abundance along the Sacramento River.They also had areas of high abundance in San Pablo Bay.

Long term changes
While B. longirostris and Acanthocyclops spp.have experienced overall declines in abundance over time, A. sinensis has mostly increased, although over a shorter time period.
The declines of B. longirostris and Acanthocyclops correspond to noted regime shifts and overall plankton declines across many species (Winder and Jassby 2011).The increase of A.
sinensis could be related to expansion following its recent invasion as it fills niches left by declining species.
The change to the seasonal pattern of B. longirostris abundance was puzzling.Prior to 1990, Bosmina experienced two peaks, one in the spring and a second peak in the fall.
Something happened in the late 1980s to cause a crash in B. longirostris abundance during the fall.This fall peak may have been an important food source for juvenile Delta Smelt.
While we may never be able to tell exactly why this peak disappeared, we discuss three theories here: changes to water management, invasive clams, and invasive copepods.
Operation of the State Water Project and Central Valley project can significantly change flows in the San Joaquin River and the channels of the South Delta, where Bosmina is most abundant (Jassby et al. 2002;Jassby 2005).Exports from these water projects cause a decrease in residence time in the South Delta, particularly during the fall peak of B.
longirostris (Hammock et al. 2019).Decreased residence time limits phytoplankton production, as well as directly exporting phytoplankton and zooplankton (Jassby et al. 2002;Hammock et al. 2019).However, this pattern first became apparent during increases to exports in the late 1970s (Hammock et al. 2019), well before the disappearance of the fall peak of B. longirostris, so the export explanation is unlikely to be the main factor driving the decrease, though it could be part of the story.
Invasive species may be a more likely explanation for the change in the seasonal peaks of B. longirostris.The entire zooplankton community experienced a decline in abundance in the western portions of the estuary due, in part, to the introduction of the overbite clam (Potamocorbula amurensis) in 1986 (Kimmerer and Orsi 1996).The overbite clam has an extremely high filtration rate, greater than any native benthic grazer, and caused a crash in chlorophyll-a in Suisun Marsh and Suisun Bay.Zooplankton suffered from loss of food resources as well as directly through grazing on zooplankton nauplii (Kimmerer and Thompson 2014;Kimmerer and Lougee 2015).Grazing rates of overbite clams peak in summer and fall as biomass increases and higher water temperatures cause increases in metabolic rates (Crauder et al. 2016), so they would be well positioned to reduce the B.
longirostris fall peak without impacting the spring peak.Could this have caused a reduction in fall B. longirostris?We feel this is unlikely to be the main cause, because B. longirostris chiefly occurs in the freshwater reaches of the Delta, whereas the overbite clam occurs mainly in areas with salinity > 2ppt (Crauder et al. 2016).While clams may not be the obvious answer, there was another important species introduction in the late 1980s with greater potential to directly impact B. longirostris.The calanoid copepod P. forbesi was introduced in 1987 and quickly became the most abundant calanoid in the system (Orsi and Walter 1991).P. forbesi peaks in abundance from July through October, overlapping with the historic peak in B. longirostris abundance, and they occur in high abundances in the South Delta where B. longirostris is common (Kayfetz and Kimmerer 2017).P. forbesi may be competing with B. longirostris for food resources during the fall when previously most other zooplankton peaked earlier in the year.
15 and low salinities of 0.15.The also exhibited a unique winter-spring abundance peak in the two highest salinities that disappeared by 2005.They invaded in 1993 (Orsi and Ohtsuka 1999), thus this pattern could reflect them settling into their ecological niche over time.
dextrilobatus typically peaks in summer (Bollens et al. 2014); however, when first introduced it also peaked in March and April (Hooff and Bollens 2004).Potential competition from T.
dextrilobatus could thus have contributed to adult A. sinensis losing its spring peak in those higher salinities.
Interestingly, A. sinensis juveniles had increasingly wide seasonal abundance peaks over time, driven in part by differing timing of the two highest salinity bins.Abundance peaks were regularly 1-2 months later in the highest salinity (16.575) than in any of the lower salinities, which all peaked around the same time.The abundance peak of the highest salinity also grew relatively larger compared to the lower salinities over time, which led to an overall widening of the seasonal abundance peak for A. sinensis juveniles.However, the width of the abundance peak in each lower salinity level also seemed to widen over time.This demonstrates shifting phenology of A. sinensis, which could be caused by changes in the timing and location of reproduction, predation, or feeding.The zooplankton community has undergone many shifts over the history of this dataset (Orsi and Ohtsuka 1999;Winder and Jassby 2011), with A. sinensis potentially having its own impacts on lower trophic level zooplankton (Kayfetz and Kimmerer 2017).Since A. sinensis is predatory it could be following the abundance shifts of other species, resulting in changes to the location and timing of reproduction, thus impacting both adult and juvenile abundances.Alternatively, this shifted phenology could be caused by changes to the salinity field in recent years which have resulted in higher salinity intrusions further into the estuary due to increased and persisting droughts (Ghalambor et al. 2021).

Conclusions
Many of the fishes in the estuary rely on zooplankton for at least part of their life cycle.Changes in prey resources can affect higher trophic levels by reducing the amount of available food or shifting the timing of peak abundance, thereby creating a mismatch between critical fish life stages and their prey.We found long-term shifts in all three of our study species.These shifts included changes in seasonality, salinity preference, and long-term abundance.Further studies investigating these patterns in additional species would be important to understand the past dynamics of zooplankton in the estuary.These results increase our understanding of the zooplankton community, which could inform the development of food web models and be matched to trends in fish abundance to examine the direct influence of declining zooplankton species on managed species.

Figure 2 .
Figure 2. Map of the study area, depicting all sampling stations.Survey abbreviations are as follows: EMP = Environmental Monitoring Program, 20mm = 20-mm Survey, FMWT = Fall Midwater Trawl, STN = Summer Townet, FRP = Fish Restoration Program.The EMP survey has some non-fixed entrapment zone (EZ) stations that move with the salinity field and are depicted here with increased transparency.

Figure 3 .
Figure 3. Seasonal patterns in Bosmina longirostris adult abundance with 99% credible intervals.Each plot represents predictions from different position on the smoothed yearly trend and may not represent exact conditions for that particular year.Predictions for different salinity values are represented by line and shading color, which is on the log scale.The y-axis limits differ among plots to facilitate comparison of seasonal trends.Absolute trends in abundance over time are represented in Fig. 4. Missing values (e.g., the months of Jan, Feb, and Dec in 1975) represent gaps in the raw data.

Figure 4 .
Figure 4. Yearly patterns in Bosmina longirostris adult abundance with 99% credible intervals.Each plot represents the pattern for a separate month.Predictions for different salinity values are represented by line and shading color, which is on the log scale.The y-axis limits differ among plots to facilitate comparison of longterm trends.Seasonal trends in abundance are represented in Fig. 3. Missing values represent gaps in the raw data.

Figure 5 .
Figure 5.Estimated values of the clustered station random intercepts for Bosmina longirostris adults.Stations within 1 km were clustered into groups.Point color indicates whether each station cluster has higher or lower B.

Figure 6 .
Figure 6.Seasonal patterns in Acanthocyclops spp.adult abundance with 99% credible intervals.See Fig. 3 for a full description.

Figure 7 .
Figure 7. Yearly patterns in Acanthocyclops spp.adult abundance with 99% credible intervals.See Fig. 4 for a full description.

Figure 8 .
Figure 8.Estimated values of the clustered station random intercepts for Acanthocyclops spp.adults.See Fig. 5

Figure 9 .
Figure 9. Seasonal patterns in Acartiella sinensis adult abundance with 99% credible intervals.See Fig. 3 for a full description.

Figure 10 .
Figure 10.Yearly patterns in Acartiella sinensis adult abundance with 99% credible intervals.See Fig. 4 for a full description.

Figure 13 .
Figure 13.Yearly patterns in Acartiella sinensis juvenile abundance with 99% credible intervals.See Fig. 4 for a full description.

Figure 14 .
Figure 14.Estimated values of the clustered station random intercepts for Acartiella sinensis juveniles.See Fig. 5 for a full description.
, but adults were first counted in samples in 1994 and juveniles first counted in 2006.Thus, adult data were filtered to the start date of 1994 and juveniles to the start date of 2006.For the other species, only adult data were available at the taxonomic resolution of our analysis.The final dataset had 33,255 samples for B. longirostris adults, 32,026 samples for Acanthocyclops spp.adults, 17,189 samples for A. sinensis adults, and 11,224 samples for A. sinensis juveniles.