Community Patterns and Environmental Associations for the Early Life Stages of Fishes in a Highly Transformed Estuary

Estuaries are important fish nursery areas, yet little is known of how environmental forcing influences estuarine fishes during their early life stages. We analyzed environmental and larval-juvenile fish community data in the upper San Francisco Estuary (SFE) from spring to early summer 1995–2017, to better understand drivers of spatiotemporal community patterns in this highly modified estuary. We evaluated community patterns based on the relative abundance and diversity of native and introduced fish in the SFE and their predominant distribution (pelagic, demersal). The upper SFE experienced a downward trend of freshwater outflow and upward trends of temperature and salinity intrusion. Fish relative abundance only showed long-term downtrends for native and introduced pelagic fish groups. The most influential habitat components for relative abundance and diversity of fish groups were in decreasing order: temperature, salinity, Secchi depth, bottom depth, and zooplankton biomass. Early life stages of native and introduced fishes were generally segregated spatially and temporally, with native fishes more associated with cooler, saltier, and higher turbidity habitats during early to mid-spring compared to introduced fishes during late spring to early summer. Community ordination showed that environmental (temperature, salinity, outflow, Secchi depth, and zooplankton biomass) and spatiotemporal factors (month and depth), explained nearly 40% of the total variance. Our results suggest that the shorter duration of planktonic and nektonic stages of demersal fish groups results in higher resiliency compared to pelagic fishes. The declining abundance of pelagic fishes overall seems to be linked to drought effects and human-induced synergistic interactions intensified by climate change.


Introduction
Estuaries are one of the most productive ecosystems on the planet and provide various ecological services such as fish nursery habitats (Barbier et al. 2011;Whitfield 2019), but they experience multiple anthropogenic impacts, including human-mediated species introductions (Cohen and Carlton 1998), hydrological and morphological modifications, habitat loss, pollution, and hypoxia (Emmett et al. 2000;Lotze et al. 2006;Hughes et al. 2015). Recent estuarine ecosystem shifts have been directly linked to human activities (e.g., Winterwerp and Wang 2013;Chevillot et al. 2017). Such shifts can be abrupt and lead to step population changes (e.g., Cloern 2007;IEP 2015), synergies among extinction drivers (Brook et al. 2008), and altered community structure and interactions (Castillo 2019). Hence, there is a need to understand ecosystem trajectories or states to inform management options (Townend 2004;Moyle et al. 2010;Mahardja et al. 2021). Fish communities have been increasingly used to help evaluate the ecological status of estuaries (e.g., Yañez-Arancibia 1985; Whitfield 2019) and can aid in the development of environmental and ecological quality objectives and standards to promote desired ecological outcomes, such as conservation of native fishes (e.g., Whitfield and Elliott 2002;Townend 2004;Barbier et al. 2011). Fish in estuaries have been grouped in several ways, with classifications generally based on particular traits within the fish community, such as life-history strategy (e.g., Elliott et al. 2007;Teichert et al. 2017), feeding and reproductive guilds (e.g., Elliott et al. 2007), salinity zones (e.g., Strydom et al. 2003), and fish morphology (e.g., Aguilar-Medrano et al. 2019). Given the large number of fish introductions in some estuaries (Emmett et al. 2000;Matern et al. 2002), and the likely influence of climate change on species invasions and interactions , the abundance and diversity of native and introduced species could provide valuable insights to interpret community patterns for non-coevolved assemblages (Moyle et al. 1986). Moreover, since the behavior of fishes during their transition from plankton to nekton influences their dispersal and retention (Leis 2006), there is a need to understand how community patterns for early life stages of fishes differ between species that remain in the pelagic zone and those who move to the demersal zone at later life stages.
The San Francisco Estuary (hereafter SFE; Fig. 1) has been profoundly transformed since the mid-nineteenth century (Nichols et al. 1986;Whipple et al. 2012) and is emblematic of an estuary facing competing ecological and socioeconomic demands. Changes reported in the upper SFE since the late twentieth century include: decreased water quality and reductions of freshwater inflow by more than half (e.g., Nichols et al. 1986;Cloern 2019); increases in the number and abundance of introduced animals and plants (Moyle et al. 1986;Matern et al. 2002;Khanna et al. 2018); contaminants (Brooks et al. 2012) and toxic cyanobacterial blooms (Lehman et al. 2013). Although multispecies declines of fish populations in the upper SFE have been reported since the late twentieth century (e.g., Bennett and Moyle 1996;Mahardja et al. 2021;Tempel et al. 2021), the populations of pelagic species reached unprecedented low levels since the early 2000s, including native and introduced fishes, and invertebrates (Cloern 2007;Mac Nally et al. 2010). Several factors were implicated in this decline, including habitat loss and degradation, water diversions, and entrainment losses at all trophic levels (Jassby et al. 2002;IEP 2015;Hammock et al. 2019). In addition, introduced species such as the clam Potamocorbula amurensis have greatly altered pelagic and 1 3 benthic food webs (e.g., Mac Nally et al. 2010;Slater and Baxter 2014). Species declines were further exacerbated by droughts (Castillo et al. 2018;Mahardja et al. 2021). For example, native pelagic fishes in the SFE like the endemic delta smelt (Hypomesus transpacificus) and longfin smelt (Spirinchus thaleichthys) moved closer to extirpation during the 2012-2016 California drought, despite both already being listed as threatened or endangered species under state and federal laws (e.g., Hobbs et al. 2017;Tempel et al. 2021). Thus, an increased understanding of the drivers of community dynamics may help to bridge the knowledge gaps between species conservation and ecosystem management.
Community analyses for early life stages of fishes of the upper SFE have been conducted in shallow bays (Meng and Matern 2001), leveed channels and sloughs (Feyrer and Healey 2003), littoral habitat (Grimaldo et al. 2004), and marshes (Matern et al. 2002;Whitley and Bollens 2014). Early life studies for fish species in open waters of the upper SFE have revealed how freshwater flow influenced the geographic distribution of species (Dege and Brown 2004), with temperature and turbidity predicting the occurrence of some pelagic fishes (Mahardja et al. 2017;Peterson and Barajas 2018). Yet, no community studies have examined the abundance and diversity for larval and juvenile stages of fishes across open waters of the upper SFE.
Evaluating the extent to which reported fish declines reflect overall trends for the early life stages of fishes at the community level is a pressing need. How the abundance and diversity of fishes in the upper SFE change in response to environmental variation is of particular importance. Yet, given the heterogeneous species composition and the large differences in abundance for early life stages of fishes in this estuary (e.g., Dege and Brown 2004;Mahardja et al. 2017), further analyses are needed to interpret community patterns. Considering the role of estuaries as major fish nursery habitats (Yañez-Arancibia 1985; Whitfield 2019), and the critical role of early life stages on the year-class strength of fishes (Houde 2008), these patterns may have major implications for restoration strategies, and climate adaptation (e.g., Quiñones and Moyle 2013;Knowles et al. 2018).
We examined the fish community in open waters of the upper SFE using larval and juvenile fish monitoring data from 1995 to 2017 based on a long-term data set (IEP 2021; Fig. 1). The goals of the study were to determine (a) the interannual environmental patterns during the period in which early life stages of fishes are present; (b) the extent to which fish abundance and diversity vary across areas, years, and habitat components; and (c) fish community patterns along environmental and spatiotemporal gradients. We show that improved interpretation of fish community patterns is possible by accounting for the predominant habitat zone of native and introduced species throughout their entire life history (pelagic, demersal), and by integrating species patterns at different levels of resolution (assemblage, dominant species, and fish groups). We then reveal emergent community patterns and environmental drivers that may be used to inform the management and conservation of native fishes in the upper SFE.

Study Area
The SFE is the largest estuary in the west coast of the USA and is the interface between ocean waters entering San Francisco Bay and the Sacramento and San Joaquin rivers, the major sources of freshwater of the SFE (Nichols et al. 1986). The upper SFE is comprised of a series of embayments and a tidally influenced network of interconnected channels which includes leveed agricultural islands over an area of 2,800 km 2 (hereafter the Delta; Fig. 1). Downstream of the confluence of the Sacramento and San Joaquin Rivers, the upper SFE includes shallow bays (Suisun Bay, Honker Bay, Grizzly Bay) surrounded to the north by large contiguous marsh (Suisun Marsh), the Napa River, and San Pablo Bay (Fig. 1). Freshwater inflow is the main driver of estuarine dynamics in the upper SFE, including transport processes, stratification and mixing, salinity, sediment and nutrient input, and dilution of pointsource pollutants (Kimmerer 2004;Brooks et al. 2012). Freshwater inflow varies greatly, both intra-and interannually, and is influenced by climate oscillations  and water diversions (Cloern and Jassby 2012). On a regional scale, the mean annual inflow into the upper SFE (ca. 1080 m 3 s −1 ) corresponds to less than 0.5% of the inflow into other large US west coast estuaries such as the Columbia River and Puget Sound (Emmett et al. 2000). Compared to other California estuaries, the SFE has a large tidal area with a salinity range of 0-25 (Monaco et al. 1992). Delta outflow is the net source of freshwater to the SFE after subtracting water diversions and depletions within the Delta and is represented by the net flow at Chipps Island (Cloern and Jassby 2012) (Fig. 1). The distance (km) along the main axis of the estuary from the mouth of the estuary (Golden Gate Bridge) to the upstream location where the near-bottom tidally averaged salinity reaches 2 is locally known as "X2," and is an ecologically useful indicator of freshwater outflow in the upper SFE (Jassby et al. 1995). During periods of high outflow, X2 is typically located in Suisun Bay (ca. 55 to 65 km), and during low outflow periods X2 can be upstream of the river's confluence (> 78 km), but its historical upstream limit has generally been near 100 km (Jassby et al. 1995).

Data Collection and Handling
Larval and juvenile fishes were collected by the California Department of Fish and Wildlife from the upper SFE every other week, from spring to early summer 1995-2017 (https:// wildl ife. ca. gov/ Conse rvati on/ Delta/ 20mm-Survey). This long-term monitoring includes up to nine estuary-wide sampling events (or "surveys") per year (IEP 2021;Tempel et al. 2021). We analyzed monitoring data at 51 stations grouped into 10 areas throughout the upper SFE (Fig. 1). Forty of these stations were sampled routinely since 1995, six stations in the north Delta were sampled since 2008, and five stations in San Pablo Bay were only sampled during high flow years (1995, 2006, 2011, and 2017). All available data from the nine surveys were considered. A total of eight surveys were conducted in the years 1995, 1996, 1999, and 2001-2004. Each survey across the upper SFE was typically completed in 1 week.
The sampling gear consisted of two joint nets, one to collect larval and juvenile fish and one to collect zooplankton. The fish net was mounted to a weighted, rigid frame set on skis, and was 4.93 m long with a mouth area of 1.51 m 2 , and constructed of 1600-μm knotless nylon mesh. One flowmeter (General Oceanics-GO, Model 2030R) was mounted across the mouth of the fish net to estimate the towing distance (m) and the volume (m 3 ) of water sampled. Zooplankton samples were collected using a Clarke-Bumpus (CB) net mounted to the top of the fish net frame. At each station, the sampling gear was lowered to near the base of the water column and three 10 min stepped-oblique fish tows were conducted. Zooplankton and water-quality samples were collected at each station, generally concurrently with the first fish tow. The zooplankton net targeted meso-zooplankton (e.g., copepods and cladocerans, and included a second GO flowmeter mounted across the mouth of the net to estimate the volume sampled. The CB net had a 160-μm knotless nylon mesh, measured 78 cm long with a mouth area of 101 cm 2 . All fish collected were identified to the lowest possible taxonomic level, with 55 taxa identified, including 54 species and one genus (Table S1). Zooplankton samples were subsampled to identify at least 200 organisms per sample (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003) or about 6% of each sample (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). Water-quality variables recorded at each station included temperature (°C), electrical conductivity (S m −1 , YSI Model 30; temperature adjusted; converted to salinity), Secchi depth (cm), and turbidity (NTU, Hach Model 2100P). Hydrological variables included Delta outflow (m 3 s −1 ) and X2 (km) during March to July 1995-2017 and were estimated from DayFlow, a two-dimensional model of historical Delta hydrology computed by California Department of Water Resources (https:// data. cnra. ca. gov/ datas et/ dayfl ow).
Based on fish swimming speed vs. fish size (He 1993;Froese et al. 2011), and the towing speed of the fish net (mean ± SD; 1.03 ± 0.1 m s −1 ), we estimate that fish larger than 122 mm are less efficiently captured by the fish net. Hence, only individuals up to 122 mm fork length (FL, if forked) or total length (TL, if not forked) were included in the analyses (Table S1). Adult fish caught in the surveys were excluded based on the size at first maturity (Table S2).
We grouped individual sampling events into "sample units," where each unit contained data averaged across stations for a given area, survey, and year. This enabled the improved representation of non-dominant species (Mahardja et al. 2017), and regional descriptions of environmental and community patterns (Castillo et al. 2018) compared to station data alone. A total of 2000 sample units were available during the period of 1995 to 2017. Fish relative abundance was expressed as catch per unit volume (CPUV; n per 10,000 m 3 ). Zooplankton was quantified both in terms of CPUV (n per m 3 ) and biomass per unit volume (BPUV, mg C per m 3 ).

Interannual Environmental Patterns
We evaluated seven environmental variables to determine potential unidirectional trends in the upper SFE over the period early life stages are present (Delta outflow, salinity, X2, temperature, turbidity, zooplankton CPUV and BPUV). Trends were evaluated using a modified Kendall test to account for seasonal covariance (Hirsch and Slack 1984) and were cross-validated using randomized tests (proportion of 100 runs with p-values ≤ to those in observed trends) (after McCune and Mefford 2018).

Fish Abundance and Diversity Across Areas, Years, and Habitat Components
Distribution of fish assemblage across water-quality variables -We compared differences in the use of three abiotic habitat components (salinity, turbidity, and temperature) between common native and introduced fishes comprising > 0.01% of the total CPUV; hence improving the representation of species distribution across habitat components in the fish assemblage (Table S1). We then used the equality of proportions test (z-test) to evaluate relative differences in the numbers of native and introduced fish taxa across three equal ranges of each waterquality variable.
Distribution, diversity, and abundance of fish groups-Native and introduced species in the upper SFE were grouped according to their predominant habitat zone throughout their life cycle (pelagic or demersal), resulting in four fish groups (introduced demersal; introduced pelagic; native demersal; native pelagic) (Table S1). Fish CPUV was log-transformed as log 10 (CPUV + 1) to reduce the influence of the most abundant species. We used Friedman's test 1 3 (T r ) to evaluate whether fish groups differed in CPUV and diversity (number of species), with areas considered randomized blocks. For each fish group, we evaluated unidirectional trends in diversity and relative abundance (summed log 10 (CPUV+1)) across years based on the Mann-Kendall test (M-K two-sided test ;Mann 1945;Kendall 1975). Trends were cross-validated as indicated for environmental variables and statistical tests were performed in SYSTAT 13 (Systat Software, Inc.).
We used nonparametric multiplicative regression (NPMR) to evaluate the extent to which the CPUV and diversity of fish groups could be predicted by five habitat components (bottom depth, salinity, Secchi depth, temperature, and zooplankton BPUV). The NPMR is a distributionfree method for fitting quantitative or binary responses, with a multiplicative kernel smoother enabling all predictors to interact. Hence, the response to any explanatory factor in the habitat model is conditioned by the values of other factors, and if any factor is intolerable, then the response is zero rather than additive (McCune 2006).
Quantitative NPMR models were run in HyperNiche (version 2.3) using the local mean Gaussian form, with the smoothing parameter (tolerance) being the span of a predictor covered by one standard deviation of the Gaussian function. The tolerance and kernel function control the portion of a multidimensional habitat space near a target point (neighborhood size). Controls on overfitting included minimum average neighborhood size (0.05 × n samples), minimum data/predictor ratio of 20, and improvement criterion for additional variables (3% for cross-validated R 2 , xR 2 ) (McCune and Mefford 2009). Randomization tests using 100 runs were conducted to evaluate whether the fit of a selected model is no better than could be obtained by chance alone, given an equal number of predictor variables and P (proportion of randomized runs with fit ≥ observed fit). To facilitate comparisons across predictors in a given model, tolerances were standardized by dividing them by the range of each predictor. Response surface plots were generated to compare patterns of CPUV and diversity across the most consistent habitat components of the fish groups. Surface response plots were based on the default minimum neighborhood size (1 for high continuity; McCune and Mefford 2009).

Fish Patterns Along Environmental and Spatiotemporal Gradients
We used canonical correspondence analysis (CCA; Ter Braak 1986) to interpret the variance in abundance of the dominant-most representative-fish species in the community, as explained by environmental and spatiotemporal covariates across the three main ordination axes. Dominant species included those taxa present in at least 20% of the sample units, and each sampling unit included at least three fish taxa (n = 1471 sample units). In addition to the covariates included in the NPMR models, we considered outflow, spatial covariates (latitude and longitude, both in decimal degrees), and temporal covariates (month and year). We used the linear combination (LC) scores of rows in the main matrix, as LC scores represent the best fit of species abundance to the environmental data (McCune and Grace 2002). Computations were performed in PC-ORD version 7.08 (McCune and Mefford 2018). Given the large differences in CPUV among dominant species, we relativized fish abundance (log 10 [CPUV + 1]) using the information function of ubiquity to emphasize the information content of species in samples and used general relativization for covariates (McCune and Grace 2002). The CCA row and column scores were standardized by centering and normalizing to unit variance so that the resulting plots were biplots (Palmer 1993). A biplot cutoff of r 2 ≥ 0.075 was used to plot covariates in the first and second ordination axes. We used Monte Carlo randomization tests (1000 runs) to test the significance of species environmental correlation.
Associations among species were evaluated in two steps using PC-ORD. First, a 2 × 2 contingency table of presenceabsence is computed for each pair of species. Species abundance is treated as presence-absence data based on the observed frequencies of "present," and "absent" cases. Expected frequencies are calculated from the marginal totals. The matrix of standardized chi-square distances (phi) is then computed among all pairs of species: where X 2 is the chi-squared statistic and N is the total observed frequency (McCune and Mefford 2018). Second, species scores (points) in the first and second axes of the joint plot are connected by lines forming a plexus diagram, with species associations set to phi > 0.5 for strong associations and phi > 0.4 and < 0.5 for weak associations (McCune and Mefford 2018).

Interannual Environmental Patterns
Spring to early summer outflows showed an interannual downtrend, with high outflows occurring only in about 20% of the years (1995, 1998Fig. 2). In contrast, X2 and salinity showed increasing peaks over time (Fig. 2), but only X2 showed an uptrend (Fig. 2). Spring to early summer average temperature also showed an uptrend over time, reaching > 19 °C during 2014-15 (Fig. 2). Turbidities were generally lower since the year 2001 but no unidirectional trend was suggested over the study period (Fig. 2). Zooplankton CPUV and BPUV varied greatly (Fig. 2) and were strongly associated (r 2 = 0.66; 21 df; p < 0.01) but neither showed unidirectional trends (Fig. 2).

Fish Abundance and Diversity Across Areas, Years, and Habitat Components
Distribution of the fish assemblage across water-quality variables-The fish assemblage showed distinct multimodal patterns in CPUV for 31 species across turbidity, salinity, and temperature (Fig. 3). In comparison to introduced species, the proportions of native species across three ranges of water-quality variables were (1) higher at low temperatures; (2) higher at high salinities; and (3) lower at low turbidities (Table 1). Distribution, diversity, and abundance of fish groups-The diversity of native fishes was generally higher in the Napa River and other areas downstream of the confluence (Fig. 4). In contrast, the diversity of introduced fishes increased at those locations upstream of the confluence. Average diversity among the four fish groups did not differ greatly in the upper SFE (F r = 5.6; 3 df; p = 0.16). The summed CPUV of native demersal fishes was lower compared to other fish groups in the upper SFE (F r = 9.24; 3 df; p < 0.05), while that of native pelagic fishes generally declined towards deeper upstream areas (Fig. 4). In contrast to demersal fish groups, pelagic fish groups showed downtrends in CPUV over the study period (Fig. 5). None of the four fish groups showed unidirectional trends in diversity across years in the upper SFE (Fig. 5).
Habitat models of CPUV and diversity for fish groups included up to four of the five covariates examined using NPMR. The most frequent habitat components across models were, in decreasing order: temperature, salinity, Secchi depth, bottom depth, and zooplankton BPUV ( Table 2). The tolerances and standardized tolerances across habitat components for fish abundance and diversity generally had comparable data bearing among models, indicating similar width of surrounding sample space (Table 2). Salinity and temperature were the most consistent habitat components to evaluate response surfaces for CPUV and diversity among the four fish groups (Table 2; Figs. 6 and 7).
Unlike introduced fish groups, the CPUV and diversity of both native fish groups were higher below 18 °C (Figs. 6 and 7). Similar to both introduced fish groups, the CPUV of the pelagic native fish group was highest at salinities below 12, while the native demersal fish group had the highest CPUV at ca. 20 (Fig. 6). Yet, for both native fish groups, the highest diversity occurred at intermediate salinities (ca. 4-12; Fig. 7a, b). In contrast, the highest diversity of both introduced fish groups occurred at salinities < 4 (Fig. 7c, d).
The only other pair of habitat components reflecting a large portion of the surface responses were Secchi depth and bottom depth, with such responses being evident only for the CPUV of native fish groups (Table 2; Fig. 8). Interestingly, both native pelagic and demersal fish groups showed highest CPUV towards the lower range of Secchi depths and bottom depths (Fig. 8).

Fish Patterns Along Environmental and Spatiotemporal Gradients
Environmental and spatial gradients reflected in the CCA ordinations generally characterized the study areas, ranging from higher salinity locations (San Pablo Bay, Napa River, West Suisun Bay) to upstream channelized areas that were associated to deeper bottom depths (e.g., south Delta, north Delta, San Joaquin River Delta), and warmer waters (e.g., south Delta; Fig. 9a). The CCA ordinations indicated the community structure for early life stages of fishes was mostly explained by five environmental factors (temperature, salinity, outflow, Secchi depth, and zooplankton BPUV) and two spatiotemporal factors (month and depth; Fig. 9b). The cumulative variance explained by the main three ordination axes was nearly 40%, with over half of such variance accounted for in the first ordination axis (Fig. 9b). Except for the native delta smelt and the introduced yellowfin goby (Acanthogobius flavimanus), native and introduced fishes were positioned along opposite sides of the 1st ordination axis, suggesting segregation between native and introduced species, with warmer months more associated to introduced fishes and cooler, earlier, months more associated to native fishes (Fig. 9b). As revealed by plexus lines, only two species were strongly associated, both introduced (threadfin shad Dorosoma petenense and American shad Alosa sapidissima). Weak associations were suggested for the remaining fishes, except for two species whose positions along the 2nd ordination axis coincided with the lowest salinities (prickly sculpin Cottus asper) and highest salinities (northern anchovy Engraulis mordax) (Fig. 9b).

Discussion
Our study on the early life stages of fishes in the upper SFE showed complementary and consistent spatiotemporal patterns at the assemblage, fish group, and dominant species levels. The observed patterns of fish distribution abundance and diversity reflect the responses of a variety of life histories, with most (57%) of the 55 fish taxa collected associated with two or more habitats (estuarine, freshwater, and marine). Although some fish species can exhibit different life histories in one or more biogeographic regions (Moyle 2002;Hobbs et al. 2019) and include a variety of migration 1 3 modes (e.g., Elliott et al. 2007;Tamario et al. 2019), only about 33% of all fish taxa collected in the upper SFE seem to have migrating life histories, including anadromy, amphidromy, catadromy, oceanodromy, potadromy, and semianadromy (Table S1). Hence, a migrating strategy may not be necessarily required by fish to exploit different habitats in transitional environments such as the upper SFE. For example, resident fish species in the North East Atlantic coast are adapted to high spatiotemporal variability of estuarine conditions (Teichert et al. 2017). Our study provides a baseline to interpret changes in fish community abundance, diversity and their environmental drivers, all which may contribute to inform recovery and conservation of species and their habitats (see "Management Implications").

Interannual Environmental Patterns
Significant long-term environmental trends were observed over the study period: decreasing outflow, and increasing X2 and temperature. These trends are consistent with the escalating severity of droughts and rising temperature partly driven by climate change (Griffin and Anchukaitis 2014; Jeffries et al. 2016). In contrast to the early warm years of our study period (1996)(1997), the more recent warmest years (2014)(2015) coincided with the highest salinities over the study period (Figs. 2 and S1), demonstrating the major influence of the 2012-2016 California drought on key habitat components for the early life stages of fishes. These findings also suggest that climate change impacts in the upper SFE (e.g., Knowles et al. 2018) may exacerbate synergistic interactions of stressors implicated in the long-term declines of native pelagic fishes in the upper SFE (Bennett and Moyle 1996;Castillo et al. 2018) (see "Influence of Habitat Changes on the Fish Community").

Fish Abundance and Diversity Across Areas, Years, and Habitat Components
Our assemblage analyses revealed distributional differences between non-coevolved species across temperature, salinity, and turbidity. These three water-quality factors were also the most influential habitat components for CPUV and diversity of fish groups. Interestingly, in the Swartkops Estuary, South Africa, these same three factors were more important in determining fish larval composition than finescale habitat factors, such as substrate type (Kisten et al. 2020). This highlights the value of those water-quality factors when interpreting fish community patterns in estuaries. Our assemblage analyses pointed out community patterns not detected for fish groups and dominant species. For example, while the increased occurrence of delta smelt and longfin smelt at intermediate and high turbidities have been demonstrated (Mahardja et al. 2017), the assemblage analyses showed that several other native species of management concern are also associated with intermediate or high turbidities, including Chinook salmon (Oncorhynchus tshawytscha), splittail (Pogonichthys macrolepidotus), and Pacific herring (Clupea pallasi) (Fig. 3).
While fish diversity did not show interannual trends for any of the four fish groups in the upper SFE, the decline in CPUV for the pelagic fish groups reflects the long-term population collapses of native pelagic fishes. Such patterns indicate the high sensitivity of pelagic fish groups to habitat degradation in the upper SFE. Although the native pelagic fish group generally showed the highest average CPUV throughout most of the study period, the introduced gobies (Tridentiger spp.) represented one of the most dominant fish taxa over the study period (Fig. 3). The three species of this genus occur in the SFE (T. bifasciatus, T. barbatus, and T. trigonocephalus), but the latter is found mostly at high salinities (Wang 2010), and is unlikely to occur in the upper SFE.
The salinity field in the upper SFE was historically more dynamic and more favorable for native species before the Delta began to be managed as a conveyance system for freshwater diversions in the late 1950s (e.g., Moyle et al. 2010;Kimmerer et al. 2019), yet the concurrent downtrends for both native and introduced pelagic fish groups suggests that degradation of the pelagic habitat has negatively impacted Fig. 3 Mean (± SD) larval-juvenile fish log-transformed catch per unit volume (CPUV; catch per 10 4 m 3 + 1) by species (bars), plotted along with the mean (± SD) turbidity, salinity, and temperature at which each species was observed in the study dataset (triangles). Species are ordered from left to right based on increasing mean value for each water quality variable. Introduced species are denoted in red font, while native species are in black ◂ Table 1 Number and percent of native and introduced larval-juvenile fishes whose weighted average catch corresponded to each range of water-quality variables in the upper San Francisco Estuary. Denoted is the z-test for the difference in proportions between native and introduced species (asterisks denote significant differences, p < 0. native and introduced pelagic species alike. In contrast, neither native nor introduced demersal fish groups display downward trends. The shorter planktonic or nektonic life stages of demersal fishes may limit their exposure to the pelagic habitat (e.g., Wellington and Victor 1989). Unlike pelagic fishes in this study that remain largely in the pelagic habitat throughout their lifespan, the transition from pelagic to benthic habitat occurs for the demersal yellowfin goby at 13 mm TL (Miyazaki 1940) and at 18 mm TL for shimofuri goby (Dotu 1958), thus potentially allowing these and other demersal species to be more resilient to environmental stressors associated with the pelagic habitat. However, it is important to note that long-lived, late-maturing anadromous native demersal species like the green sturgeon (Acipenser medirostris) and white sturgeon (A. transmontanus) have experienced long-term declines and are, respectively, endangered and of special concern due to habitat degradation and losses within and outside the upper SFE (Miller et al. 2020).
The majority of fish sampled by the survey in this study are in their larval or early juvenile life stages (Table S1); therefore, net avoidance by pelagic and demersal species at their planktonic and nektonic stages should be limited. Such potential effects are further minimized by the maximum fish length considered in our study and by excluding individuals larger than the size at first maturity (Table S2). Moreover, separate data from depth logger tests for the fish net used in this study confirmed sampling occurs over most of the water column (Trishelle Tempel, unpublished data). These field tests showed the percentage of the water column sampled (mean ± SD) was 96 ± 9% for tows at depths ≤ 10.7 m; 94 ± 8% for tows at depths ≤ 13.4 m, and 69 ± 11% for tows at depths ≥ 13.4 m. The highest depth range in those tests accounted for only 1% of the sample units in our study, with the deepest area (the Confluence) having an average depth of ca. 12 m (Fig. 4).
As in the case of juvenile-adult pelagic fishes across the upper SFE (Castillo et al. 2018), our study showed that larvaljuvenile stages of introduced pelagic fishes are mostly freshwater species and native pelagic species tend to be associated with higher salinities. We further showed that the Fig. 4 Diversity and relative abundance, summed logtransformed (n per 10 4 m 3 + 1) for four fish groups in sample units (mean ± SE) across areas. Also shown is the bottom depth (mean ± SE) per area and the distance from the mouth of the San Francisco Estuary to each area (ordered from downstream to upstream) different salinity distribution patterns between native and introduced species apply similarly to pelagic and demersal species (Fig. 7). The association between native species and higher salinity levels is likely partly driven by the increasingly unsuitable freshwater and low salinity habitats of the SFE for these species (e.g., Moyle et al. 2010;Hobbs et al. 2017). Nevertheless, based on the fish distribution across habitat components in the assemblage ( Fig. 3; Table 1) and in fish groups (Figs. 6 and 7), segregation between noncoevolved species seemed primarily due to temperature, with the highest CPUV of introduced fishes observed at the warmest conditions. Our study further showed that segregation between native and introduced fish groups across temperature and salinity was more clearly defined in terms of diversity than CPUV. This pattern could be due to diversity gradients reflecting broader changes in overall species tolerance to temperature and salinity compared to the sharper changes in CPUV, as most of the CPUV per fish group could be contributed by just a few species. This interpretation is consistent with our fish assemblage analyses showing different proportions between native and introduced species across temperature and salinity (Table 1).
Zooplankton BPUV in the upper SFE was relatively high during the 2012-2016 California drought (Figs. 2 and S2), coinciding with the precipitous declines of age-0 delta and longfin smelts (Hobbs et al. 2017). This suggests that any potential trophic benefits due to increased zooplankton may have been overridden by generally unfavorable abiotic conditions for native fishes during that period, including unusually warmer, saltier, and clearer waters (Figs. 2, S1, S2), and such conditions likely limited their co-occurrence with prey and thus their survival. Zooplankton BPUV was the least predictive habitat component among fish groups, being only associated with CPUV of introduced pelagic fishes (Table 2), with a seemingly stronger relationship at salinities below 8 (Fig. S3). The generally lower explanatory power of zooplankton BPUV in our habitat models could be due to the concurrent response of zooplankton and fish to abiotic variation, the low resolution of total zooplankton BPUV as an index of potential food supply, or both. Given that different prey items often vary in their nutritional value and predator avoidance behaviors, more specific food-based metrics (i.e., zooplankton species and their developmental stage) perhaps should be considered in future habitat models for select fish species. Since monitoring of depleted fish populations such as the native osmerids of the upper SFE (i.e., delta smelt and longfin smelt) can be difficult, the use of proxy species with a similar prey base may be useful to evaluate the combined Fig. 5 Relative abundance (summed log-transformed CPUV + 1) and total diversity (total number of species) per year for four fish groups from spring to early summer 1995 to 2017. Significant trends for Kendall τ are denoted by an asterisk (2-sided M-K test) potential of food and outflow management actions during spring months to increase the prey base important to larval and juvenile pelagic fishes.

Fish Patterns Along Environmental and Spatiotemporal Gradients
Consistent with the previous assemblage and fish group analyses, our ordination for dominant fish species supports temperature as the primary driver of segregation between native and introduced species in the upper SFE (Fig. 9). Ordination analyses further revealed that temperature was associated with month, with months reflecting differences in the timing of spawning and recruitment into the gear, where introduced species generally predominated in later and warmer months than native species. This suggests that overall fish recruitment into the survey gear may shift earlier with climate change. Additionally, outflow influenced the distribution of species along the salinity gradient in our ordination analyses. Although sampling units were generally aggregated by area, some of them showed high variability. Such highly dynamic areas are mostly found in the downstream portion of the upper SFE and seem the combined result of high hydrological variation and water management operations in the Delta (e.g., Kimmerer 2004).
Compared to native fishes, introduced fishes were generally more associated with higher values of Secchi and bottom depths, and zooplankton BPUV. The observed temporal segregation between native and introduced fishes in the upper SFE is consistent with the reported community patterns for early life stages of fishes in the south Delta (Feyrer 2004). Our analyses also suggest that such segregation is Table 2 Habitat models for relative abundance (summed log 10 (CPUV + 1)) and diversity (number of species) of four fish groups (NP: native pelagic, ND: native demersal, IP: introduced pelagic, ID: introduced demersal) as a function of habitat components across sample units.
Shown are nonparametric multiplicative regression's tolerances, standardized tolerances, average neighborhood size (N*), and cross-validated R 2 (xR 2 ; P < 0.01) a Surface plot for salinity and temperature (Fig. 6); b Surface plot for Secchi and bottom depth (Fig. 8); c Slice plot for salinity at selected zooplankton CPUV (Fig. S3); d Surface plot for salinity and temperature (Fig. 7); e Slice plot for salinity at selected Secchi or bottom depths (Fig. S4

Influence of Habitat Changes on the Fish Community
Broad water-quality monitoring across the upper SFE began in the 1970s and has revealed (1) decreases in suspended particulate matter, chlorophyll-a, and dissolved oxygen (DO); (2) increases in Secchi depth, ammonium, and nitrate + nitrite (Cloern 2019), and (3) higher water temperatures during most months of the year (Bashevkin et al. 2022). The rate at which habitat components continue changing has critical implications for future regime shifts and species extinction risks, particularly for stressors enhanced by climate change (Brook et al. 2008;Tempel et al. 2021) and factors associated with long-term habitat alteration (Nichols et al. 1986;IEP 2015). Hence, there is a need to interpret the variation in habitat components at several temporal and spatial scales. For example, although our trend analyses did not reveal a decline in turbidity in the upper SFE overall, such decline was suggested in the north Delta, south Delta, and lower San Joaquin River Delta (Fig. S2), in line with the reduction of total suspended solids in the Delta between 1998 and 2010 (Hestir et al. 2013). Consistent with the increasing CPUV of both native fish groups towards lower Secchi depths, our fish assemblage analyses showed that all but one of the native fish species predominated at turbidities > 35 NTU, underscoring the value of this habitat component for native fishes. The long-term increase of Secchi  (Knowles et al. 2018), in combination with eutrophication, algal blooms, and low outflow could further limit DO in various spots within the SFE and potentially lead to hypoxic events (Cloern et al. 2020) such as those reported in other estuaries (e.g., Hughes et al. 2015;Scanes et al. 2020). Zooplankton BPUV and CPUV showed no unidirectional trends over the study period (Fig. 2), yet important zooplankton prey for some pelagic fishes in the upper SFE seems in short supply (Slater and Baxter 2014). Zooplankton BPUV declined greatly in the early 1980s and its species composition has and will likely continue to shift . Moreover, the temporal co-occurrence of early life stages of fishes and their prey could be altered by their ontogenetic sensitivity to warmer seasonal climate (e.g., Komoroske et al. 2014;Jeffries et al. 2016). Abundance peaks for phytoplankton and some species of zooplankton in the upper SFE seem to be occurring earlier (Merz et al. 2016). This could alter the temporal overlap between juvenile fish and their zooplankton prey, as reported in the Gironde Estuary, southwest France (Chevillot et al. 2017). Given the marked differences in phenology between non-coevolved pelagic fishes in the upper SFE, and their concurrent CPUV declines, climate-driven impacts on these fish groups could include altered trophic interactions, ecophysiological limitations, and increased extinction risk (e.g., Penn and Deutsch 2022).
Our study suggests that several of the abiotic changes implicated in the upper SFE pelagic organism decline that occurred in the early 2000s (e.g., Cloern 2007;Mac Nally et al. 2010) still persist and their ecosystem impacts may continue to be intensified by climate change and other stressors. These abiotic changes likely facilitate novel  Fig. 1, and b early life stages for the 12 most common fishes in relation to covariates and ordination statistics. Introduced fishes are denoted in red font and strong and weak plexus associations among fishes are respectively denoted by blue and red lines species interactions at multiple trophic levels in the upper SFE Khanna et al. 2018;Castillo 2019), especially as new species continue to be introduced (Cohen and Carlton 1998;Mahardja et al. 2020). Projected increases in temperatures and frequency of droughts could further expand blooms of toxic cyanobacteria (Lehman et al. 2013;Ustin et al. 2014). The extent to which these changes could influence the fish community composition and interactions also depends on the ability of each species to adapt, for example by shifting their range or migrating, as reported for fishes in South African estuaries (e.g., Whitfield 2019).
Considering the arid to semiarid climate of the SFE watershed (Faunt and Sneed 2015), the projected severity of droughts in the SFE under climate change (Griffin and Anchukaitis 2014;Knowles et al. 2018) could magnify the effects of hydrological alteration on fish communities (e.g., Sakaris 2013). For instance, drought-induced salinity intrusions will likely be exacerbated by water diversions (e.g., Moyle et al. 2010;Castillo et al. 2018) and projected sea level rise (Knowles et al. 2018). Hence, the geographical distribution of fish communities and their spawning and rearing areas in the upper SFE could largely differ from their historic distributions (e.g., Moyle 2002;Wang 2010). For native fish species, particularly those dependent on the low salinity zone (0.5-6.0; IEP 2015), those changes could restrict their habitats toward the narrower and deeper upstream channels in the upper SFE, potentially resulting in lower overall fish abundance and species diversity (Figs. 8 and S4), as well as lower carrying capacity (e.g., Castillo 2019).
There may be a couple of reasons why only the introduced demersal fish group seemed to expand its upstream range during the 2012-2016 California drought (Fig. S5). First, the introduced demersal fish group has a higher percent of migratory species (57%) potentially facilitating distributional shifts compared to other groups: introduced pelagic (17%), native demersal (29%), and native pelagic (40%) (Table S1). Second, the temperature-salinity patterns across years and areas (Fig. S1) and the response surfaces across salinity and temperature (Figs. 6, 7) suggested that distributional shifts were less likely or apparent in other fish groups due to physiological limitations. The native demersal fish group was more associated with higher salinities but cooler waters than other fish groups, potentially limiting its upstream shift during droughts. Likewise, the association of the native pelagic fish group to cooler waters may have prevented an upstream shift due to physiologically inadequate warmer upstream areas despite more suitable salinities (e.g., Jeffries et al. 2016). The apparent ability of the introduced demersal fish group to tolerate and move across a broad range of salinity and temperature indicates that the species within this group may be more capable of facing the stressors associated with climate change.
Based on our study, climate change may also favor introduced pelagic fishes over native pelagic fishes in mostly warm and freshwater Delta. However, such climate effects may be less apparent compared to more immediate anthropogenic threats to native fishes (Quiñones and Moyle 2013). Yet, despite the relative dominance of introduced pelagic fishes in the Delta (Fig. 4), their overall paucity in more saline regions of the upper SFE suggests that they may not necessarily benefit from climate change-induced salinity intrusion (via sea level rise and more frequent droughts). However, saltier waters downstream of the Delta may become less suitable to native species as those areas warm (e.g., Jeffries et al. 2016), increasing the potential for new introduced species colonizing those areas (e.g., Cohen and Carton 1998). When all areas are combined over the study period, both introduced pelagic and demersal fishes had higher CPUV compared to both native fish groups at least since 2015 (Fig. 5). Though prediction on the likely future of diversity and CPUV for the four fish groups is outside the scope of our study, CPUV and diversity projections for introduced vs. native, and demersal vs. pelagic could be developed through our study findings. For example, habitat associations from studies such as ours can be integrated with downscaled climate projections (e.g., Knowles et al. 2018;Stern et al. 2020) and the projected effect of water quality changes on species physiology (e.g., Komoroske et al. 2014;Jeffries et al. 2016).

Management Implications
The observed trends for pelagic fish groups in our study are consistent with the downward trends of several species of pelagic fish in the upper SFE (e.g., Castillo et al. 2018;Mahardja et al. 2021;Tempel et al. 2021). Sustainable management of these public trust resources may require implementing long-term actions consistent with population recovery plans for listed and vulnerable species (e.g., USFWS 1996) along with more immediate near-term management actions to limit the risk of extinction (Hobbs et al. 2017). Borgelt et al. (2022) found that many data-deficient species worldwide are likely facing high extinction risk. This highlights the need to evaluate the population status of declining native species not presently classified as threatened or endangered under state and federal laws. Considering the observed and projected climate change impacts (Knowles et al. 2018), and the deteriorating habitat conditions within the rearing period of fishes in the upper SFE, both proactive and innovative management actions seem warranted to address the declining long-term resilience of pelagic fish communities in open waters (Mahardja et al. 2021) and to maintain life history diversity of critically endangered species (Hobbs et al. 2019). Effective recovery goals would likely involve large-scale restoration targeting natural abiotic habitat variability (Nobriga et al. 2005;Moyle et al. 2010;Castillo et al. 2018), including adaptive flow augmentations across multiple areas and seasons (e.g., Castillo 2019; Sommer et al. 2020); sediment supply (SFEI-ASC 2016), and reducing stressors such as contaminants (Brooks et al. 2012) and high temperatures (Ustin et al. 2014).