Peer Reviewed Title: Benthic Assemblage Variability in the Upper San Francisco Estuary: A 27-Year Retrospective Journal Issue:

We used multivariate methods to explore changes in benthic assemblage structure over 27 years (1977–2003) at four monitoring stations located along a salinity gradient in the upper San Francisco Estuary. Changes in benthic assemblage composition were assessed relative to hydrologic variability and to the presence of the high-impact invader Corbula amurensis in the estuary. We also explored the composition of benthic assemblages during a recent collapse of several pelagic populations in the upper estuary. Our results show that the Corbula invasion had both direct and indirect effects on the benthos in the estuary, causing significant changes in assemblage structure. We found no unprecedented patterns of benthic assemblage composition during the period of the Pelagic Organism Decline (2000–2003) in the upper estuary. Hydrologic variability was associated with significant changes in benthic assemblage composition at all locations. Benthic assemblage composition was more sensitive to mean annual salinity than other local physical conditions. That is, benthic assemblages were not geographically static, but shifted with salinity, moving down-estuary in years with high delta outflow, and up-estuary during years with low delta outflow, without strong fidelity to physical habitat attributes such as substrate composition or location in embayment vs. channel habitat. Organism abundance and species richness showed a bi-modal distribution along the salinity gradient, with lowest abundance and richness in the 5 to 8 psu range. We conclude that the continuity of benthic assemblages and community metrics along the salinity gradient is a powerful and necessary context for understanding historical variability in assemblage composition at geographically static monitoring stations.


Local Identifier: jmie_sfews_11014
Abstract: We used multivariate methods to explore changes in benthic assemblage structure over 27 years  at four monitoring stations located along a salinity gradient in the upper San Francisco Estuary.Changes in benthic assemblage composition were assessed relative to hydrologic variability and to the presence of the high-impact invader Corbula amurensis in the estuary.We also explored the composition of benthic assemblages during a recent collapse of several pelagic populations in the upper estuary.Our results show that the Corbula invasion had both direct and indirect effects on the benthos in the estuary, causing significant changes in assemblage structure.We found no unprecedented patterns of benthic assemblage composition during the period of the Pelagic Organism Decline (2000)(2001)(2002)(2003) in the upper estuary.Hydrologic variability was associated with significant changes in benthic assemblage composition at all locations.Benthic assemblage composition was more sensitive to mean annual salinity than other local physical conditions.That is, benthic assemblages were not geographically static, but shifted with salinity, moving down-estuary in years with high delta outflow, and up-estuary during years with low delta outflow, without strong fidelity to physical habitat attributes such as substrate composition or location in embayment vs. channel habitat.Organism abundance

INTRODUCTION
A substantial portion of the water used by California's agricultural and urban consumers is drawn from the Sacramento-San Joaquin Delta, in the upper San Francisco Estuary.As demand for this water outpaces the supply, many options for managing freshwater resources and environmental quality are being considered for the upper estuary (Lund and others 2008).Environmental and biological data, developed through decades-long monitoring, Benthic Assemblage Variability in the Upper San Francisco Estuary: A 27-year Retrospective san francisco estuary & watershed science 2 provide a unique opportunity to explore ecosystem responses to a broad range of environmental variability, and analyses of these data can provide a model for potential future plans based on what has come before.Benthic invertebrates have been an important component of the estuarine ecosystem (Alpine and Cloern 1992;Kimmerer and Orsi 1996;Nichols 1985;Orsi and Mecum 1996;Jassby and others 2002;Feyrer and others 2003), and a better understanding of the historic benthic assemblage responses to hydrologic and biological variability in the estuary could improve our ability to forecast the outcomes of proposed management and restoration activities in the upper estuary.
Previous work in the San Francisco Estuary has shown that interannual variability in freshwater flow, due to both anthropogenic and natural factors, influences variability in the abundance or survival of several biotic communities in the estuary ranging from phytoplankton to fishes (Jassby and others 1995;Kimmerer 2002aKimmerer , 2002b)).Historical relationships between flow and benthic invertebrate assemblage composition in this estuary have been demonstrated (Nichols 1979(Nichols , 1985;;Nichols and Thompson 1985;Markmann 1986;Nichols and others 1990;Alpine and Cloern 1992;Hymanson and others 1994).However, recent investigations of pelagic populations (Kimmerer 2002a;Sommer and others 2007) showed precipitous declines in abundance in several pelagic fish populations during recent years with high freshwater outflow when population increases were expected.These declines in pelagic organism populations have encouraged a series of assessments of biotic processes and pollutant effects in the estuary.It is therefore timely to re-assess benthic invertebrate assemblage patterns and their relationship to potential forcings.
The "remarkable invasion" of the Asian bivalve Corbula amurensis in the San Francisco Estuary (Carlton and others 1990; hereafter referred to as the Corbula invasion) has been associated with a cascade of ecological changes that have affected much of the tidal estuary.Vigorous filter-feeding by the clam, which has become both abundant and wide-spread in the estuary (Carlton and others 1990;Hymanson 1991), has been linked to a sharp decline in phyto-plankton biomass and a removal of seasonal phytoplankton blooms in the northern estuary (Alpine and Cloern 1992;Jassby and others 2002;Jassby 2006).Corbula amurensis is tolerant of a broad range of salinities, from 2 to 32 psu in a laboratory setting (Nicolini and Penry 2000), and occupies a broad range of substrates (Carlton andothers 1990, Hymanson 1991).It became numerically dominant in much of the upper estuary, including the reach from San Pablo Bay to the Lower Sacramento River since shortly after its discovery in the estuary in 1986.
The declines in phytoplankton have not been limited to the regions of the estuary where the clam is abundant (Jassby 2006).The suppression of chlorophyll a (a proxy for phytoplankton biomass) was first observed in the Suisun Bay region where dense populations of C. amurensis have the capacity to filter a volume equivalent to the overlying water nearly daily (Alpine and Cloern 1992;Werner and Hollibaugh 1993).The zone of depleted chlorophyll has extended well beyond the geographic distribution of the clam (Jassby 2006), as chlorophyll produced in Suisun Bay has been depleted by grazing clams, and thus is no longer available to be delivered to upstream locations through tidal sloshing (Kimmerer and Orsi 1996;Jassby and others 2002).The decline in phytoplankton biomass in the upper estuary, including the lower Sacramento and San Joaquin rivers (Jassby 2006) has led to food limitation (Muller Solger and others 2002) which has been linked to declines in zooplankton (Kimmerer and others 1994;Orsi and Mecum 1996;Kimmerer and Orsi 1996) and several pelagic-feeding fish species, including one endangered species, the delta smelt, Hypomesus transpacificus (Bennett and Moyle 1996;Kimmerer and others 2000;Kimmerer 2002a, Feyrer andothers 2003).Effects of the Corbula invasion on the benthic invertebrate assemblages in Grizzly Bay have been explored previously (Peterson 2002), but the nature of benthic assemblage changes across the salinity gradient of the upper estuary, both upstream and downstream of Grizzly Bay, has not been explored until now.
The Corbula invasion occurred at the beginning of a prolonged drought in the San Francisco Estuary watershed (1987)(1988)(1989)(1990)(1991)(1992).Since the hydrology in the estuary has subsequently returned to more average may 2010

Study Area
The upper San Francisco Estuary (Figure 1) comprises a series of embayments and channels which lead from the Pacific Ocean to an inland river delta, where the Sacramento and San Joaquin rivers join.The Sacramento River contributes roughly 80% of the freshwater input to the estuary.The San Joaquin River, which historically contributed a large volume to the estuary, is presently mostly diverted for storage and delivery to agriculture, and the remaining flows are managed for fish passage and environmental concerns.Freshwater inputs to the delta from the upstream watershed are managed by reservoir releases and water diversions.The primary goal of the management of freshwater inflow to the estuary is to guarantee that the inland delta remains flushed with high quality fresh water for export from the southern delta (Kimmerer 2002b).These management activities affect flow patterns through the estuary (Kimmerer 2002b;Jassby 2005).
Even with considerable, active management of the region's water resources, natural climate variability, including the interannual variability of precipitation and snow melt in the San Francisco Estuary's 140,000 km 2 watershed, is a dominant influence on hydrologic conditions in the estuary (Knowles 2000).Total inflow can vary by a factor of 10 among years and interannual variations in inflow lead to variations in transport of many aquatic constituents as well as variations in the location of the low-salinity zone in the estuary (Jassby and others 1995;Kimmerer 2002b;Kimmerer 2004).Seasonal to interannual variability in freshwater inflow and the resulting variability in the estuarine salinity field is an important factor in structuring the benthic assemblages of the estuary (Nichols 1979;Nichols and Thompson 1985;Peterson 2002).

Environmental Data
Monitoring of benthic invertebrates and water quality has been conducted in the upper San Francisco Estuary, including the Sacramento-San Joaquin Delta, under the auspices of the California conditions it has become clear that historic relationships between hydrology and population densities of key pelagic fish species have failed to predict the severe population declines observed in recent years.Severe declines in fish populations have been observed since 2000 (Bennett 2005;Feyrer and others 2007;Sommer and others 2007).The decline of pelagic fish populations in the upper San Francisco Estuary has been referred to as the "Pelagic Organism Decline" (POD) (Sommer and others 2007).Hypotheses posited to explain the changing relationships between fish populations and hydrologic conditions in the estuary have emphasized food web changes attributed to invasive species (primarily C. amurensis) but also have implicated increased occurrences of toxic Microcystis algae, changes in pesticide use, and alterations in water management practices which could be hampering the recovery of fish populations (Sommer and others 2007).
Nearly three decades of monitoring of the aquatic environment and benthic invertebrates of the upper San Francisco Estuary as part of the Interagency Ecological Program's monitoring effort (http://www.water.ca.gov/iep/activities/emp.cfm, Environmental Monitoring Program) has chronicled changes in assemblage structure during extreme drought and flood conditions in the estuary, as well as during the establishment of several invasive species.Because monitoring was conducted at several locations along the salinity gradient of the upper estuary, we can compare responses in benthic assemblages from multiple salinity regimes and habitat types to physical and biological perturbations.
In this paper we examine the underlying model of benthic community variability in the context of environmental forcings across multiple sites.We explore the nature of benthic assemblage variability in the upper estuary.We also examine how the structure of the benthic assemblages has been affected locally by interannual variability in hydrology as well as the invasion of Corbula amurensis, (even at locations where C. amurensis is not found).Finally, we explore differences in benthic assemblage structure before and during the POD for indications of changes in the benthic assemblage that might help identify conditions or processes that are leading to declines in pelagic fish populations.
were then used to compute annual flow statistics (Figure 2A).Hydrologic variability in the watershed of the San Francisco Estuary results in a broad range of conditions in the estuary's aquatic habitats.Years with prolonged high outflows result in lower than average salinities at most upper estuary locations as well as transport, dilution and mixing that are associated with sustained elevated flows.Years with either short pulses of high flow followed by prolonged low outflow or simply prolonged low outflow result in elevated salinities in much of the upper estuary and a different flow, mixing and transport regime.In order to assess the effects of hydrologic variability on the benthos we sought to identify hydrologic year types and assess benthic assemblage response to the conditions presented in wettest and driest of years.For the purposes of this study, the distribution of annual Department of Water Resources (CDWR), in cooperation with the Interagency Ecological Program (IEP) since the mid 1970s.IEP data from several program components were used for our analyses.
Hydrological data were obtained from the Dayflow database administered by CDWR (http://www.water.ca.gov/dayflow/).The delta outflow (Dayflow parameter "QOUT") was selected to represent general hydrologic conditions in the upper estuary for the purposes of this study."QOUT", the "net delta outflow at Chipps Island" (Figure 1), is a computed estimate, produced by CDWR, of net flow of fresh water from the Sacramento-San Joaquin Delta, seaward past Chipps Island.It is derived by performing a water balance computation based on net inflows of fresh water to the delta and tidal forcing in the upper estuary.Daily flow data (obtained on March 20, 2007) were used to compute monthly means which IEP salinity data, reported as specific conductance in µmhos/cm, were converted to practical salinity units (psu) using the 6th-order polynomial equation given in Schemel (2001), based on Lewis (1980).Annual statistics were computed from mean monthly values.
Sediment data were provided by the IEP Environmental Monitoring Program (Karen Gehrts, CDWR, unpublished data).Particle-size analysis and determination of the percentage of organic matter (primarily peat soil and vegetable detritus) were conducted on one-liter soil sub-samples collected during benthic invertebrate sampling events.Sediment analyses were conducted by DWR Soils Laboratory, Bryte, CA.Grain size and organic matter analyses are reported separately (CDWR unpublished data).Annual statistics for each constituent were computed from monthly values.
The BEST function in PRIMER was used to rank correlations between the benthic assemblage composition and environmental data.Environmental variables were normalized prior to constructing a dis-similarity matrix using Euclidean distance.The dissimilarity matrix based on environmental variables was compared to a Bray-Curtis dissimilarity matrix computed from fourth-root transformed annual average organism abundance data, as described below.The abiotic variable which best matched the pattern of benthic assemblage composition was salinity (correlation = 0.77), the next best was for salinity and sediment percent fines (correlation = 0.73), and the addition of further variables resulted in reduced correlation values.Thus, the analyses that follow focus on patterns of assemblage composition with regard to salinity and percent fines.

Benthic Assemblage Data
Of 21 benthic monitoring stations sampled by the IEP monitoring program, four stations situated within a broad salinity gradient and located in four major regions of the upper San Francisco Estuary have yielded the longest and most consistent records: D41A-C in San Pablo Bay (SPB), D7-C in Grizzly Bay (GRB), D4-L in the lower Sacramento River (LSR), and D28A-L in Old River (OLR) (Figure 1).San Pablo Bay, where station SPB is located, is a shallow, polyhaline embayment and represents the seaward, most saline end of the salinity gradient explored in this study.Grizzly Bay where station GRB is located is a shallow sub-embayment of the mesohaline Suisun Bay complex.The lower Sacramento River station (LSR), which is generally oligohaline habitat, is located at the mouth of the Sacramento River, near its confluence with the San Joaquin River and at the western edge of the inland Sacramento-San Joaquin Delta.The Old River station (OLR) lies along a natural channel of the San Joaquin river drainage which is currently one of the main pathways through which fresh water is conveyed from the Sacramento River to the southern delta to be exported from the estuary to agriculture and urban users in other regions of the state.
IEP benthic sampling was conducted seasonally to monthly from 1977 through fall of 2003.At each of the four benthic monitoring stations, three or four grab samples per sampling event were taken from a boat using a Ponar dredge which samples an area of 0.053 m2.Samples were then washed through a Standard No. 30 stainless steel mesh screen (0.595 mm openings).All material remaining on the screen after washing was preserved in approximately 20% buffered formaldehyde containing Rose Bengal dye.In the laboratory, benthic invertebrates were sorted from each sample, identified to lowest possible taxon (species where possible), counted, and preserved in 70% ethanol.Identification and enumeration of invertebrates was conducted by Hydrozoology, Newcastle, California, for the entire period of record.Details for the monitoring program and program data can be found at http://www.baydelta.water.ca.gov/emp/Metadata/benthos_metadata.html.
The period of record for each of the stations discussed in these analyses spans nearly 30 years.Some variation in benthic sampling frequency did occur among stations and over time.For stations GRB, LSR, and OLR, spring and fall sampling began in 1977, with monthly sampling occurring from June 1980 to October 2003.Station SPB was originally part of a separate program (Schemel and others 1988(Schemel and others , 1989(Schemel and others , 1995)), and was sampled bi-monthly from 1987 to mid-1991.Beginning in September 1991, when it became part of the IEP monitoring program, SPB was sampled at the same monthly interval as the other stations (GRB, LSR, and OLR).The availability of benthic assemblage data is shown in Appendix A. After October 2003 sampling at all four stations was cut to quarterly, and previous analysis (Peterson and Vayssières 2006) showed that the change in sampling frequency produced a systematic under-sampling of several species that confounded analyses of assemblage changes beyond 2003, and thus the later data were excluded form this study.
Mean monthly averages of benthic species abundance were computed using data from the samples collected each month at each station.Abundance data from the 2 to 12 sampling events per year were aggregated to compute annual mean species assemblage composition for each station-year combination.Prior to analyses, Balanus improvisus (which is considered an encrusting species rather than soft-bottom benthos), and unidentified species were removed.Data from dominant oligochaete species Limnodrilus hoffmeisteri and Varichaetadrilus angustipenis were aggregated to eliminate the effects of changes in taxonomic identifications over time on the analyses, as the two abundant species were distinguished from each other only after the mid-1980s.
Multivariate analyses were carried out on benthic invertebrate assemblage data from the four sampling sites using PRIMER software (Clarke and Gorely 2006).Assemblages similarities were explored using non-metric multi-dimensional scaling (NMDS; Clarke and Green 1988;Clarke 1993) on Bray-Curtis dissimilarity indices computed from 4th-root transformed annual mean invertebrate abundance data.The Kruskal stress formula (Kruskal 1964) was computed to show goodness of fit of two-dimensional data plots to multi-dimensional dissimilarity data, and only plots with stress ≤0.1, indicating good to very good fit, were used in this study.Cluster analysis, using group average clustering, was also performed on the Bray-Curtis dissimilarities; to further assess the relationships between the benthic assemblages found in our samples.
The early portion of the record (1977)(1978)(1979) was valuable because it included the last year of an extreme drought.Because there were as few as two samples per year in this portion of the record, we had to consider limiting our analyses to two samples per year (May or June and September or October) for the entire study which would give even coverage for all years, versus using all available data (up to 12 sampling events per year) which would maximize the information we could gather from assemblage patterns for the 24 years with nearly monthly data.We compared the two Bray-Curtis dissimilarity matrices resulting from these two approaches; using results from 999 randomly permuted samples (Clarke and Warwick 1994; the "RELATE" function in PRIMER-6 software; Clarke and Gorely 2006) to test the hypothesis of no relation between Bray-Curtis dissimilarity matrices from the two sets of data.The dissimilarity matrices are used in this test because they are at the core of subsequent analyses.The analysis resulted in the rejection of the null hypothesis that the two matrices were unrelated (ρ = 0.98, significance p = 0.01; of 999 permutations run, zero permutations may 2010 resulted in ρ ≥ 0.98).Because the rank correlations showed that the two dissimilarity matrices were reasonably similar, subsequent analyses of benthic community change were completed using all available data.
Community metrics such as species richness and organism density were computed from the biotic data to examine trends in community metrics along the salinity gradient and over time.Community metrics that were computed for years where less than five samples were collected have been omitted from the trend analyses.Annual data computed from years with less than five samples yielded slightly inflated (approximately 10% greater) annual organism density and greatly reduced (approximately 20% less) species richness, and these values weren't suitable for comparison with the remaining data.
NMDS plots provide a visualization of the amongsample Bray-Curtis dissimilarity (dissimilarity equals 1-Bray-Curtis similarity) by placing samples that are most similar close together and less similar samples at an increasing distance.The computation of Bray-Curtis dissimilarity takes into account the presence and abundance of all species in each sample, in this case, the 4th root transformed mean abundance of each of the species found at a site during a year.While the distribution of samples (represented by symbols) in a NMDS plot implies a pattern of benthic assemblage composition, the symbols generally do not convey the species composition of the samples being compared.We have enhanced the NMDS plots by substituting pie-diagram representations of benthic assemblage composition for the symbols that were initially plotted.We used the X and Y coordinates of the symbols from NMDS plots to spatially organize the pie diagrams, which display untransformed annual average abundance of the numerically dominant benthic species for each sample.The size of the pie diagrams is scaled to represent the average of monthly abundance of benthic invertebrates for each station-year sample.
An analysis of similarity (ANOSIM) procedure (Clarke and Green 1988;Clarke 1993) is a multivariate analog of ANOVA which compares the average of rank similarities among replicates within a group to average of rank similarities arising from all pairs of replicates between groups.By randomly reassigning the groups to which replicates belong and recalculating the rank similarities, then comparing the true rank similarity (R, reported as ρ) to those found from random permutations, a probability of randomly achieving the result can be calculated and thus a significance level (p) assigned to the rejection of the hypothesis that the similarities between and within sites are equal.The ANOSIM analysis was applied to Bray-Curtis dissimilarity indices to analyze the significance of differences in faunal composition at each of the four monitoring stations between groups of samples: (1) samples in wet versus dry years; (2) samples prior to the Corbula invasion (1977-1985, preinvasion) and years since (1987-2003, post-invasion); and (3) samples taken during the Pelagic Organism Decline (POD) years [2000][2001][2002][2003] and all remaining years since the Corbula invasion (1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999).
Groups of samples that were determined to be significantly different using ANOSIM were further analyzed using similarity percentage analysis (SIMPER; Clarke 1993).The SIMPER routine assesses the contribution of each species to between-groups dissimilarity and reports average abundance, the average betweengroups dissimilarity, variability of the dissimilarity contribution among samples, and the relative contribution of each species to the between-groups dissimilarity.The SIMPER analysis was conducted on 4th-root transformed data.To summarize the effective between-groups differences in abundance for all species that contributed at least 2 percent to between-groups dissimilarity in the SIMPER analyses, the mean and standard error of untransformed species abundances were calculated from samples belonging to four categories: pre-Corbula invasion and dry years, pre-Corbula invasion and wet years, post-Corbula invasion and dry years, and post-Corbula invasion and wet years.
To determine whether significant differences in benthic invertebrate population density occurred between pre-and post-Corbula invasion assemblages, one-way Analyses of Variance (ANOVA) were performed on 4th-root transformed mean annual benthic organism density from two time periods (1977-1985 and 1987-2003).Prior to the ANOVA homogeneity of variance of the abundance data was assessed using Levene's test (Levene 1960).Where more than two groups were compared a Tukey multiple comparison (Tukey 1953) at the 0.05 significance level was used to discern significant differences among group means.Analyses were conducted in Minitab statistical software.
The monitoring stations are located along a salinity gradient.The highest salinities were found at SPB, while OLR remained "fresh" (less than 1 psu) throughout the study period (Figure 2).Salinity variability was least in dry years at SPB and greatest in wet years.Conversely, the variability of salinity was least in wet years at LSR and OLR and more pronounced in drier years.At GRB, salinity variability was high during all but the most hydrologically extreme years.
Fine sediments were common at all four stations.While both SPB and GRB were primarily composed of silt and finer sediments, LSR and OLR sediments contained variable proportions of sand and fines within and among years (Figure 3A).At LSR the mean percentage of fines tended to be highest in high delta outflow years including 1983and 1995-1998and lowest during low outflow years including 1977, 1987-1991, and 2002 (Figure 3A) (Figure 3A).The organic content of sediments at SPB and GRB was uniformly low over the study period (Figure 3B).The organic frac-tion of sediments was generally below 20% at OLR, but reached up to 50% at LSR where organic matter content was generally negatively correlated with percent fines (r (27) = -0.43,p < 0.05) over the study period.

Variations in Benthic Assemblage Structure Among the Four Locations
The organism density and species richness of benthic assemblages varied among stations as well as among years within stations (Figure 4).During the study period, organism density (abundance m -2 ) was generally greatest at SPB and LSR (grand mean organism densities are 14,800 m -2 and 11,100 m -2 , respectively) and least at OLR and GRB (grand mean organism densities are 5,000 m -2 and 3,700 m -2 , respectively).Abundances at LSR were significantly less after the Corbula invasion than they had been in the preinvasion period (F (1, 26) = 22.7, p < 0.05).Other sites do not show significant reductions in mean organism abundance during 1987-2003 compared to the 1977-1985 period.Species richness is highest on average at SPB (grand mean = 19) followed by OLR, LSR, and GRB (grand mean 16, 14, and 9 species per sampling event, respectively).
Benthic assemblages tended to be more similar within stations than among neighboring stations, but there was occasional inter-station overlap in assemblage composition, generally during hydrologic extremes.Both multi-dimensional scaling plots (Figure 5) and cluster analysis (Appendix B) of Bray-Curtis dissimilarity of among station-year samples from all of the four long-term monitoring sites from 1977 to 2003 demonstrate the predominance of within-station similarity.Information from both analytical techniques is used to illustrate most accurately the relationships between samples (Figure 5).
The NMDS analyses used here present a two-dimensional representation of multidimensional dissimilarities, in this case, benthic assemblage composition among station-years.There are some instances where station-years that have very small dissimilarity may appear to be more spatially distant from each other than to other adjacent but more dissimilar stationyears.To better visually distinguish groups of samples that had less than 50% dissimilarity in the cluster analysis (Appendix B), they are delineated in the NMDS plot (Figure 5), and referred to as "Group" 1 to Substrate composition at four long-term benthic monitoring stations, including grain size distribution (A), summarized as percent fines (particles < 0.075 mm), the remainder being sand or small gravel; and (B) percent combustible (at 400°C) organic matter, primarily peat soil and vegetable detritus.The two components are not to be considered additively.For the box plots, top and bottom "x" represent 5th and 95th percentiles, the box spans the 25th to 75th percentiles, line across the box is the median, and the small square in the box is the annual mean.
11 in the text.The goodness of fit of the two-dimensional representation in the NMDS plot to the multidimensional dissimilarity is measured by the Kruskal stress formula, which is ≤ 0.1, indicating good to very good fit, for Figure 5 and for the subsequent NMDS visualizations used in this paper.
Instances of inter-station overlap of benthic assemblage composition are evident during the study period, primarily in samples from GRB and LSR (Figure 5).During the extreme drought in 1977, species assemblages from GRB and LSR became more similar to each other than they were to other assemblages found at their respective stations (Figure 5, Group 3).In 1982, 83, and 84, during and immediately following a prolonged high-outflow period, the benthic assemblages found at GRB, were more similar to assemblages from LSR from 1978 than to other GRB samples (Figure 5, Group 4).
Assemblage structure varied along the salinity gradient in the estuary, with salt-tolerant species dominating the assemblage at SPB, and salt-intolerant species dominating in the low salinity at OLR.The average Bray-Curtis dissimilarity between Group 11, which is composed entirely of SPB samples, and Group 2 which is composed entirely of OLR samples (Figure 5) is 99%, indicating that there was nearly complete species turnover between assemblages from the two locations.Dissimilarity coefficients between neighboring stations range from 56 to 67 percent dissimilarity, based on pair-wise SIMPER comparisons of assemblage composition between stations, indicating species overlap among neighboring locations along the salinity gradient.
The distribution of benthic assemblages along the X-axis of the NMDS was strongly correlated with salinity (R 2 = 0.97) along an exponential decay line (Figure 5B).When community metrics from the benthic invertebrate samples were organized along the X-axis of the NMDS, bi-modal patterns of organism density and species richness were evident (Figure 6).Plotting values of community metrics using the NMDS X-values instead of salinity allowed us to show greater detail in species richness and organism density patterns in the low salinity habitats, whereas these patterns would be distorted and compacted in the very low salinities, below 2psu, by an X-axis based on salinity.Based on these visualizations (Figure 6) two things are readily apparent: first, the decrease in organism abundances at LSR in post-Corbula invasion years (after 1986) and, second, species richness is generally greater in post-invasion years than pre-invasion years.

San Pablo Bay (SPB)
The  7).The amphipod Corophium heteroceratum and the bivalve Mya arenaria were also commonly abundant (Figure 7).While there was no clear long-term trend in total organism abundance at SPB from 1987 to 2003, mean annual benthic invertebrate abundance tended to be least in wet years, including 1995, 1996, and 1998 (Figure 4).The highest mean annual abundance, however, occurred in 2003 which was a wet year for the purposes of this study.
An analysis of similarities (ANOSIM) of benthic invertebrate assemblage composition in wet vs. dry years at SPB showed that there was a significant difference in the assemblage structure between the two year types (ρav = 0.60, p = 0.003).Analysis of the species that dominated each assemblage type and the species that contributed most to the contrast between them (using SIMPER) showed that ostra-cod Eusarsiella zostericola, amphipods A. abdita, C. heteroceratum, Monocorophium acherusicum, and Ampelisca lobata, polychaetes Nephtys cornuta franciscana, Euchone limnicola, and Sabaco elongatus, bivalves Theora lubrica, and Musculista senhousia, and tunicate Molgula manhattensis were more abundant during dry years (Table 1A).During wet years C. amurensis and N. hinumensis tended to be more abundant than in dry years and the polychaete Marenzelleria viridis was regularly present in the assemblage (Table 1A).
A second ANOSIM indicated that the benthic assemblage composition observed in POD years was not significantly different from composition in all other post-invasion years (ρav = 0.03, p = 0.3).Because there was no significant difference between the assemblages in the two time periods no SIMPER analysis was conducted.

Grizzly Bay (GRB)
The species composition and population density of the GRB benthic invertebrate assemblage was highly variable over the 1977-2003 study period (Figure 4).During the drought year 1977, mean annual organism density was the highest observed at this station over the study period, and numerically dominant species included Mya arenaria, Ampelisca abdita and Monocorophium acherusicum, bivalves and amphipods common in the SPB assemblage (Figure 8).Mean annual organism density was also high at GRB during the 1982-1983 period of sustained high delta outflow, when taxa including Corbicula fluminea and Americorophium stimpsoni, commonly found at LSR, dominated the GRB assemblage.After the Corbula invasion (1986), there was a clear change in the structure of the benthic assemblage, with C. amurensis, the cumacean Nippoleucon hinumensis, and the amphipod Corophium alienense remaining numerically dominant in the assemblage in all post-invasion years (Figure 8).ANOSIM between invertebrate assemblages found in wet and dry years revealed a significant difference in assemblage composition between the two year types (ρav = 0.37, p = 0.003).Likewise, there was a significant difference between assemblage compositions found in pre-versus post-Corbula invasion years (ρav = 0.89, p = 0.004).SIMPER comparisons of assemblage composition in wet versus dry years prior to the Corbula invasion, showed that the benthic assemblage at GRB was dominated by amphipods Americorophium stimpsoni and A. spinicorne, bivalve Corbicula fluminea, and mermithid nematodes in wet years, and by the bivalves Mya arenaria, Macoma petalum, Musculista senhousia, amphipods A. abdita, Monocorophium acherusicum, Grandidierella japon-   1B).In dry years since the invasion, crustaceans Nippoleucon hinumensis, A. abdita, G. japonica, as well as polychaete worms Pseudopolydora kempi, Streblospio benedicti, Heteromastus filiformis, oligochaete Tubificoides brownae and several other species have generally also been abundant (Table 1B).
C. amurensis has been numerically dominant at GRB in all hydrologic conditions since the clam invaded the upper estuary in 1986, though on average, it has been more abundant at this station in dry years (Table 1B).
SIMPER comparison of the assemblage composition found that in wet years pre-versus post-Corbula invasion show that the post-invasion wet assemblage included the addition of C. amurensis and N. hinumensis, increases in the abundance of C. alienense, M. viridis, G. daiberi and the oligochaete worm Tubificoides heteroceratum, as well as reductions in the abundance of amphipods A. stimpsoni, A. spinicorne, bivalve Corbicula fluminea, and several oligochaete, nematode, and nemertean worm species (Table 1B).The comparison of dry year assemblage composition pre-versus post-invasion showed that the addition of C. amurensis and N. hinumensis, as well as a large reduction in the average abundance of bivalves M. arenaria, Musculista senhousia, Macoma petalum, amphipod M. acherusicum and several other species contributed to dissimilarity between the postinvasion dry assemblages (Table 1B).
Analysis of assemblage patterns during POD years versus other post-invasion years using ANOSIM showed that the assemblage composition during POD years was not significantly different from all other post-invasion years (ρav = -0.19,p = 0.96).

Lower Sacramento River (LSR)
The benthic invertebrate assemblage composition at LSR was variable over the study period.In most years, the amphipods Americorophium stimpsoni and A. spinicorne, the bivalve Corbicula fluminea and oligochaete worms Varichaetadrilus angustipenis and Limnodrilus hoffmeisteri were numerically dominant in the benthic assemblage (Figure 9).In 1977In and 1987In -1992, when delta outflow was low, polychaete worms Boccardiella ligerica and Laonome sp. were dominant in the assemblage (Figure 9).The highest mean organism density was observed in 1979 and 1980, lowest in 1996 and 2001 (Figure 4E).There was a clear decrease in organism density at LSR over the study period (Figures 4E and 9).Mean organism density was significantly greater from 1977-1986 than from 1987-2003 (F (1, 26) =22.7, p < 0.05).
Comparison of benthic assemblage composition between wet and dry years using ANOSIM showed that there were significant differences between species assemblages found in different hydrologic year types (ρav = 0.43, p = 0.008), likewise significant differences were found between species assemblages found before and after the Corbula invasion (ρav = 0.79, p = 0.004).In wet years prior to the Corbula invasion, the LSR benthic assemblage was dominated by the amphipods A. stimpsoni, A. spinicorne, the bivalve Corbicula fluminea, and oligochaete worms (Figure 9; Table 1C).In dry years prior to the Corbula invasion the polychaete Boccardiella ligerica was abundant and numerically dominant.Amphipods Grandidierella japonica, Monocorophium acherusicum and Melita nitida, isopod Synidotea laevidorsalis were also present and abundant in dry years (Table 1C).
During wet years after the invasion, amphipods A. stimpsoni, A. spinicorne, and Gammarus daiberi, polychaetes Laonome sp. and Marenzelleria viridis as well as oligochaete and mermithid nematode worms were present in the assemblage (Table 1C).B. ligerica and C. amurensis were regularly present during dry years after the Corbula invasion (Figure 9; Table 1C).
The observed changes in wet year type assemblages since the Corbula invasion included reductions in the abundance of A. stimpsoni, A. spinicorne, and C. fluminea and the additions of G. daiberi, Laonome sp., M. viridis, C. amurensis, N. hinumensis, and snail Melanoides tuberculata to the assemblage (Table 1C).The observed changes in dry year assemblages since the Corbula invasion include the additions of C. amu-rensis, G. daiberi, Laonome sp., N. hinumensis, a turbellarian flat worm, and the oligochaete Ilyodrilus templetoni to the assemblage, as well as reductions in the average abundance of C. fluminea (Table 1C).
Comparison of assemblage composition during POD years vs. all other post-invasion years using ANOSIM showed that there was no significant difference between assemblage composition from the two periods (ρav = 0.01, p = 0.42).

Old River (OLR)
The OLR benthic assemblage was the least variable of the four stations in terms of species composition.
There was no obvious trend in organism abundance, but there was a general increase in species richness over the study period (Figure 4).The invertebrate assemblage at OLR was dominated by the amphipod Americorophium stimpsoni, the bivalve Corbicula fluminea, oligochaete worms Varichaetadrilus angustipenis, Limnodrilus hoffmeisteri, and the polychaete Manayunkia speciosa in most years (Figure 10).In some years the amphipods Gammarus daiberi or Americorophium spinicorne were also abundant here.The highest mean annual organism density was observed in 1979 and the lowest in 1977.There was no sampling at OLR in 1978 Comparison of benthic assemblage composition between wet and dry years using ANOSIM showed that there were significant differences between species assemblages found in different hydrologic year types (ρav = 0.16, p = 0.04).There was, likewise, a significant difference in the composition of assemblages found before and after the Corbula invasion nematodes, the oligochaetes Limnodrilus hoffmeisteri and Varichaetadrilus angustipenis, polychaetes M. speciosa, and Laonome sp., ostracod Cyprideis sp., gastropod Melanoides tuberculata, and bivalve Pisidium casertanum during wet years and slightly elevated abundance of amphipods A. spinicorne, A. stimpsoni, and G. daiberi, as well as oligochaete Branchiura sowerbyi during dry years (Table 1D).
The most consistent changes in the wet year type assemblage after the Corbula invasion included the additions of the polychaete Laonome sp., amphi-(ρav = 0.73, p = 0.004).SIMPER analyses shows that the dissimilarity between wet and dry year types pre-invasion was primarily due to differences in species abundance rather than species composition (Table 1D).Dominant species, including the bivalve C. fluminea, amphipods Americorophium stimpsoni and A. spinicorne, and oligochaetes were generally more abundant in the wet year type (Table 1D).After the Corbula invasion, the most consistent contrast between assemblages found in wet versus dry year types was, again, from increased abundance of several species, including Dorylaimid and Actinoliamid The pie-charts show untransformed species abundances.Groups of years, shown within ovals are more similar (>65% Bray-Curtis similarity) to ring-mates than to nearest neighbors.Species and corresponding colors are organized by taxonomic groups; mollusks are represented in blue shades, crustaceans in gold and brown, annelids are represented by polychaetes in purples, oligochaetes in greens, and other species are in remaining colors.Wet and dry years are indicated by "W"and "D," respectively.pod G. daiberi, and gastropod M. tuberculata to the assemblage, as well as increased abundance of M. speciosa, and reduced abundance of oligochaete Limnodrilus udekemianus, as well as amphipod A. stimpsoni (Table 1D).Consistent changes in the dry year type assemblage that have been observed after the invasion include decreased abundance of oligochaete Limnodrilus udekemianus and amphipod Americorophium stimpsoni (Table 1D).
Comparison of assemblage composition during POD years vs. all other post-invasion years using ANOSIM showed that there was no significant difference between assemblage composition from the two periods (ρav = 0.06, p = 0.3).

Benthic Assemblage Response to Hydrologic Variability
The four long-term stations lie on a salinity gradient and salinity variations at each station were a function of hydrologic conditions (Figure 2), primarily delta outflow, a relationship that has been explored in detail by others (Jassby and others 1995;Monismith and others 2002;Kimmerer 2004).Variations in benthic assemblage composition were strongly associated with changes in salinity (Figure 5).Observation of benthic assemblage changes over the 27-year study period indicated that changes in benthic assemblage composition along the salinity gradient were continuous, despite discontinuities in physical habitat between stations (Figure 5).In years of hydrologic extremes assemblage compositions commonly found at each monitoring station generally followed salinity.Assemblages moved down-estuary in years with high delta outflow, and up-estuary during years with low delta outflow (Figure 5), without strong fidelity to substrate type or channel vs. embayment habitat, physical conditions which were most strongly contrasted between adjacent stations GRB and LSR (Figures 1 and 3).
With assemblage data organized by NMDS X-value (a proxy for salinity) rather than by geographic location or time, continuous patterns in community metrics were also evident (Figure 6).The bimodal distributions of abundance and species richness along the salinity gradient are similar to those presented by Remane and Schlieper (1971) whose model of species richness postulated that neither marine nor freshwater species are well adapted to the mesohaline salinity range and that only the few truly euryhaline-adapted species can thrive in the 5 to 8 psu salinity zone.Salinities at GRB are most often within that 5 to 8 psu range, and species richness is indeed generally least in GRB samples (Figures 4 and 6).During the highest outflow years (1982, 1983, 1995, 1996 and 1998) salinity at SPB reached this 5 to 8 psu range as well, and organism density is very low at this station in 1995, 1996, and 1998 (Figures 4 and 8) (Figures 4 and 8).There were no benthic data available from SPB during the high outflow years 1982 and 1983.There was one year in the study period when salinity at LSR was in the 5 to 8 psu range, 1977, and the population density in this year was not clearly below average for this station.This sample may not fit the expected pattern for several reasons, including the low number of samples for 1977, and the downward shift of mean abundance at LSR due to the significant decline in mean abundance that occurred after the Corbula invasion.Remane and Schlieper's (1971) model of a species minimum in the low-salinity zone between 5 and 8 psu was based on work in the non-tidal Baltic, and several studies have been conducted in other brackish water habitats since then which explore the pattern (Wolff 1974;Deaton and Greenberg 1986;McLusky and others 1993;Ysebaert and others 1993;Yamamuro 2000;Attrill 2002).While the bimodal distribution of organism density and species richness over the salinity gradient have been shown for many estuarine and non-tidal brackish water systems, the shape of the distribution can be affected by disturbances such as organic or chemical pollution.For instance, the rebound of species richness in the tidal fresh water of the Rhine-Meuse Estuary, the Netherlands, was dampened in the late 1960s by organic and chemical pollution, which limited the fauna of this region to a few pollution-tolerant oligochaete species (Wolff 1973).The upstream fauna was likewise limited in the Schelde and Ems estuaries, the Netherlands (Ysebaert and others 1998) and in the James River Estuary, USA (Diaz 1989).In the Forth Estuary in Scotland, McLusky (1989) reported that while species richness in the tidal freshwater remained low due primarily to organic pollution (McLusky and others 1993), organism abundance was very high in the upper reaches of the estuary as a result of organic enrichment.A decade later, however, after a reduction in pollution, the distributions of organism abundance and species richness in the upper Forth Estuary more closely resembled those shown in the present study.Likewise, an improvement in species richness has been found in the tidal fresh water of the Rhine-Meuse Delta in response to restoration (Smit and others 1997).Though it is a nutrient rich, urbanized estuary, most of the tidal fresh water of the upper San Francisco Estuary, particularly near sites LSR and OLR, does not suffer from over-production of phytoplankton or anoxia, which could inhibit a broad range of benthic invertebrates.Thus, increased abundance and richness in the tidal fresh water is apparently not dampened in the upper San Francisco Estuary by these common effects of pollution.
While the relationship between benthic assemblage composition and salinity is not a new idea, our findings highlight the importance of considering estuarine benthic assemblages as a continuum occurring over the salinity gradient when assessing changes in benthic assemblages at individual sites.Our analysis of benthic community composition arranged along a salinity continuum reveals the dynamic nature of assemblage as a function of the geographic distribution of salinity conditions in the estuary.These results demonstrate that conceptualizing benthic assemblages in the continuum fashion is a powerful and necessary approach for understanding historical variability in assemblage composition at the geographically static monitoring stations.We believe that the continuum concept could also be useful in anticipating the likely outcomes in the estuary of future hydrologic changes due to both climate change and changes in water management strategy.Additionally, these findings imply that trends in assemblage-metrics at a sampling location depend, in a continuum fashion, on salinity and thus on hydrologic conditions.Diversity measures or biotic indices used to assess trends in the "health" or "condition" of benthic assemblages or, ultimately, the local aquatic environment, often rely at least in part on community measures such as abundance and species richness (Karr and Chu 1999).The relationship between outflow, salinity, and community metrics shown here implies that such indicators need to be salinity insensitive or salinity corrected to function properly.
The ordination of the benthic assemblage data from all four monitoring stations over the study period resulted in a necessarily abstract model of benthic assemblage change over time and in various salinity regimes.A more detailed view of assemblages on a station by station basis displays the community compositions that gave rise to the broader patterns discussed above.At each of the four monitoring stations, assemblages from similar hydrologic year types tended to cluster together due to similarities in dominant organisms and overall organism density, with assemblages from wet years tending to appear toward one side of the NMDS plots and assemblages from dry years toward the other (Figures 7, 8, 9 and 10).The contrast in assemblage composition between wet and dry year types is most evident at GRB (Figure 8) where a transition from an assemblage dominated by Mya arenaria, Monocorophium acherusicum and Ampelisca abdita (all of which were generally found at SPB) in 1977, to assemblages dominated by Americorophium stimpsoni, Corbicula fluminea and oligochaetes (which were generally found at LSR) in 1982-1984 after a period of sustained high outflow, is most striking.Also evident in the station by station ordinations are the stark changes in assemblage composition that occurred after the Corbula invasion in 1986, and these are most apparent at GRB and LSR (Figures 8 and 9).

Changes Associated with the Corbula amurensis Invasion
There were clear differences between pre-and post-invasion assemblage compositions at all four monitoring locations, but each location's benthic assemblage responded somewhat uniquely to the invasion.C. amurensis became established at three of the four monitoring stations.At LSR, C. amurensis was not a dominant factor in assemblage structure.However, observed changes in the benthic assemblage at that station after the Corbula invasion were likely a response to ecosystem-level changes associated with the invasion.Species richness was generally higher in the post-invasion period at all of the stations where pre and post-invasion data are available (Figure 6).At all of the stations some of the benthic assemblage change between pre-and post-invasion periods was attributed to the increased frequency and abundance of new species, C. amurensis among them, that have invaded the estuary during the course of the study period.
Significant differences were found between the pre and post invasion assemblage compositions at all stations except San Pablo Bay where insufficient data were available for the pre-invasion period.However, supplemental data exist which describe the assemblage in years before the clam arrived.Data collected by USGS in the San Pablo Bay shallows in 1975, 76, and 86 (Janet Thompson, USGS, unpublished data) showed that before the invasion, Ampelisca abdita was numerically dominant, and the bivalves Mya arenaria, Gemma gemma, as well as polychaete Streblospio benedicti were also commonly abundant.While A. abdita remained numerically dominant after the Corbula invasion the other species were rarely abundant.C. amurensis and Nippoleucon hinumensis, along with A. abdita dominated the post-invasion assemblage.Overall organism densities reported in earlier studies were similar to those found in recent years such as 2003.
The upper San Francisco Estuary, from the Grizzly Bay region eastward, has had a clear reduction in phytoplankton biomass since 1987 (Jassby and others 2002;Jassby 2006).The loss of phytoplankton is believed to be due, in part, to the increase in benthic filtration of phytoplankton by C. amurensis after its establishment (Alpine and Cloern 1992;Jassby and others 2002).Significant downward trends in chlorophyll a (a proxy for phytoplankton biomass) have been reported by Jassby (2006) for GRB, LSR, and OLR, but no significant decline was reported for SPB (Jassby 2006).At LSR the downward trend in phytoplankton biomass is coincident with significant reductions in overall benthic organism abundance.While there was not a concurrent significant decline in organism abundance at GRB, there was a clear and persistent shift to numeric dominance by C. amurensis in the benthic assemblage there after the invasion (Figure 8).At OLR the changes in assemblage composition were primarily due to increases in abundance and the addition of species which were recent invaders to the system, and no clear change in organism density was observed that would reflect a response to a reduction in phytoplankton.
Species richness was generally higher in the postinvasion period at all stations where comparisons can be made (Figure 6).Among the highest contributing factors to assemblage dissimilarity between the assemblages was the presence and increased abundance of species that were recent invaders to the San Francisco Estuary (Table 1).Coincident with the addition of C. amurensis to GRB was the stable, numerically dominant accompaniment of several additional invasive species, including Corophium alienense, Nippoleucon hinumensis, Marenzelleria viridis and increased occurrence of Gammarus daiberi.Both C. alienense and G. daiberi had been identified in the estuary well before the Corbula invasion (in 1973 and1983, respectively;Cohen and Carlton 1995).It is unclear what factors led to the persistence and abundance of C. alienense in particular only after the mid-1980s.At LSR, the post-invasion period was associated with increased abundance of several recent invaders, including G. daiberi, Laonome sp., M. viridis, N. hinumensis, and M. tuberculata.While M. tuberculata was first identified in the estuary in 1972 (Cohen and Carlton 1995), the remaining species were likely introduced contemporaneously with or after C. amurensis; as the first records in the San Francisco Estuary for Laonome sp., M. viridis, and N. hinumensis, were 1989, 1991, and 1986respectively (Cohen and Carlton 1995).C. amurensis did not occur at OLR, though increased occurrence and abundance of G. daiberi, M. tuberculata, and Laonome sp. have been found in the assemblage regularly in the postinvasion period.
Corbula amurensis is by far the most notable recent invader due to its broad distribution and sheer abundance in the estuary, not to mention the grazing pressure it has asserted on estuarine phytoplankton.However, the clam is just one of many new species making its home in the estuary.Each of the new invasive species mentioned above is commonly found at two or more of the monitoring stations shown in this study, implying that they are well adapted to the variable euryhaline environment in the upper estuary.That new invasive species would play such a large role in the changing nature of assemblage composition is not surprising, however, because the preponderance of invasives and the press of new invaders in the San Francisco Estuary has been extraordinary (Cohen and Carlton 1995;Cohen and Carlton 1998).What is surprising is that these five recent invaders (Corophium alienense, Nippoleucon hinumensis, Marenzelleria viridis, Gammarus daiberi, and Laonome sp.) have remained fairly stable in the assemblages and that they have remained numerically dominant since the establishment of C. amurensis.

Benthic Assemblage Structure During the Recent Declines in Pelagic Fish Populations
Assemblage structure during the years of the Pelagic Organism Decline (POD) was not significantly different from other post-invasion years at any of the stations.This indicates that there was nothing unusual about the composition of benthic assemblages during the POD at the four long-term monitoring stations, such as a loss of a species or group of species, which might result from a change in the environment during that period.Further, several hypotheses commonly proposed as potential explanations for the POD were not supported by our analysis.
First, while the decline in overall benthic organism density at LSR concurrent with the Corbula invasion supports the hypothesis that suppression of phytoplankton biomass by the clam has affected populations that rely on pelagic supplies of organic carbon, as many previous studies have shown in detail (Alpine and Cloern 1992;Mueller Solger and others 2002;Jassby and others 2002;Jassby 2006), we found no evidence from the benthic abundance data that the influence of benthic grazing underwent a significant change coincident with the POD.
Next, no decline in amphipod species (which along with other crustaceans are important prey for pelagic fish) was evident during the POD.Such a decline might have indicated the possibility of toxicity effects due to increased use of pyrethroid pesticides in the watershed.It is important to note, however, that the monitoring sites used in this study are removed from urban and agricultural areas near which acutely toxic sediment and water samples have been collected (Werner and others 2000;Werner and Moran 2008).Further, because the introduction of pyrethroid pesticides was preceded by the use of organophosphorus pesticides, and because the relative toxicities of these two pesticides are similar among life forms (in other words, arthropods are most sensitive, worms are least sensitive, mollusks are also generally not very sensitive), it would be difficult to ascribe community-level effects of the change in pesticide use by measuring changes in relative abundance of the various life forms in the community.Nonetheless, we feel that the persistence of amphipods, at least at the IEP long-term monitoring sites, implies that the role of pesticides in the POD may be limited.
Finally, no decline in snail and bivalve species concurrent with the POD was detected.Such a trend could have indicated that microcystines from toxic algal blooms which have become more common in the estuary in recent years (Lehman and others 2008), had affected benthic invertebrate populations.Species-specific effects of Microcystis blooms include massive die-offs of aquatic snails (Thiaridae, Physidae, and Planorbidae) during Microcystis blooms (Kotak and others 1996;White and others 2005).Bioaccumulation of the toxic microcystine toxin was found to occur in snails (Kotak and others 1996;Zurawell andothers 1999), bivalves (Yokoyama andPark 2002), mussels (Williams 1980), clams (Prepas and others 1997) and crayfish (Lirås and others 1998).We believe the lack of a decline in these species coincident with the POD indicates that microcystines probably did not have a broad effect in the upper estuary.However, we expect that communitylevel effects are more likely to be found in waters nearer to the source of the toxin.

CONCLUSIONS
Previous studies in the San Francisco Estuary have described benthic species assemblages in a geo-may 2010 graphically static manner (Lee and others 2003), or examined changes in benthic assemblage structure as related to environmental variability at a specific locality (Nichols 1985;Nichols and Thompson 1985;Nichols and others 1990).Nichols (1985) and Nichols and others (1990) pointed out the influence of hydrologic forcing in shifting the benthic assemblage composition and community metrics.Our study, however, examines changes in benthic assemblage composition at four adjacent stations positioned along the salinity gradient over a prolonged study period that encompassed broad variety of hydrologic conditions.With this breadth and depth of data we demonstrate that benthic assemblages were not geographically static, but shifted geographically with salinity, under the influence of hydrologic forcing, without strong fidelity physical habitat attributes such as substrate composition or location in embayment vs. channel habitat.Moreover, as assemblages moved geographically, the values of community metrics such as organism abundance and species richness stayed with the assemblage on the journey, leading to a perceived change at the geographically static monitoring stations.We feel that it is necessary, therefore, to consider the concept of continuity of assemblage composition along the salinity gradient of the upper estuary when making assessments of changes in benthic assemblages over time.
Clearly, other factors affected benthic assemblage composition, and in the San Francisco Estuary, invasion was an important one.The Corbula invasion had both direct and indirect effects on the benthos in the estuary, causing significant changes in assemblage structure.Though we expected to find changes in the benthos during the period of the Pelagic Organism Decline (2000)(2001)(2002)(2003) in the upper estuary that might help us understand the environmental factors contributing to the decline, we found no unprecedented patterns of benthic assemblage composition during that time.
The fact that we found robust, continuous patterns in the benthos from four monitoring locations that are so broadly spaced is a testament to the value of consistent long-term biological and environmental monitoring.It was only by capturing environmental variability and ensuring taxonomic consistency over a multi-decadal scale that we were able to capture robust patterns of biotic change.

Figure 2
Figure 2 Hydrologic conditions and salinities at the four monitoring locations 1977-2003.Mean monthly Delta Outflow Index (A) and mean monthly salinity values for the four long-term benthic monitoring stations (C, D, E, F) are shown plotted by year in box charts, where top and bottom "x" represent 5th and 95th percentiles, the box spans the 25th to 75th percentiles, line across the box is the median, and the small square in the box is the mean of monthly flow values for the year.Hydrologic year type categories (B) are shown as "wet" (W), and "dry" (D).
Figure 3Substrate composition at four long-term benthic monitoring stations, including grain size distribution (A), summarized as percent fines (particles < 0.075 mm), the remainder being sand or small gravel; and (B) percent combustible (at 400°C) organic matter, primarily peat soil and vegetable detritus.The two components are not to be considered additively.For the box plots, top and bottom "x" represent 5th and 95th percentiles, the box spans the 25th to 75th percentiles, line across the box is the median, and the small square in the box is the annual mean.

Figure 4
Figure 4 Trends in organism density and species richness from four long-term benthic monitoring stations.The number of samples (N) are shown by year below each plot.

Figure 5 Figure 6
Figure 5 NMDS of benthic assemblages from the four longterm monitoring stations 1977-2003.The NMDS plot (A) is based on Bray-Curtis dissimilarities of assemblage composition among all samples.Numerals shown on the plot represent the year (yy) of the sample and the four stations are shown in four shapes and shades.Groups of years, shown within ovals have dissimilarities <50%.The fitted curve (B) shows relationship between mean annual station-year salinity and the X-axis value of the sample in the NMDS plot.The curve fits to an exponential decay with R 2 = 0.97.

Figure 7
Figure 7Benthic assemblage structure inSan Pablo Bay from 1987-2003.Pie charts representing dominant species abundance and relative organism density (pie diameter) of benthic assemblage taxa are arranged based on NMDS analysis of Bray-Curtis dissimilarities.The pie-charts show untransformed species abundances.Groups of years, shown within ovals are more similar (>65% Bray-Curtis similarity) to ring-mates than to nearest neighbors.Species and corresponding colors are organized by taxonomic groups; mollusks are represented in blue shades, crustaceans in gold and brown, annelids are represented by polychaetes in purples, oligochaetes in greens, and other species in remaining colors.Wet and dry years are indicated by "W"and "D," respectively.

Figure 8
Figure 8 Benthic assemblage structure in Grizzly Bay from 1977-2003.Pie charts representing dominant species abundance and relative organism density (pie diameter) of benthic assemblage taxa are arranged based on NMDS analysis of Bray-Curtis dissimilarities.The pie-charts show untransformed species abundances.Groups of years, shown within ovals are more similar (>65% Bray-Curtis similarity) to ring-mates than to nearest neighbors.Species and corresponding colors are organized by taxonomic groups; mollusks are represented in blue shades, crustaceans in gold and brown, annelids are represented by polychaetes in purples, oligochaetes in greens, and other species in remaining colors.Wet and dry years are indicated by "W"and "D," respectively.

Figure 9
Figure 9 Benthic assemblage structure in the lower Sacramento River from 1977-2003.Pie charts representing dominant species abundance and relative organism density (pie diameter) of benthic assemblage taxa are arranged based on NMDS analysis of Bray-Curtis dissimilarities.The pie-charts show untransformed species abundances.Groups of years, shown within ovals are more similar (>65% Bray-Curtis similarity) to ring-mates than to nearest neighbors.Species and corresponding colors are organized by taxonomic groups; mollusks are represented in blue shades, crustaceans in gold and brown, annelids are represented by polychaetes in purples, oligochaetes in greens, and other species in remaining colors.Wet and dry years are indicated by "W"and "D," respectively.

Figure 10
Figure 10 Benthic assemblage structure in Old River from 1977-2003.Pie charts representing dominant species abundance and relative organism density (pie diameter) of benthic assemblage taxa are arranged based on NMDS analysis of Bray-Curtis dissimilarities.The pie-charts show untransformed species abundances.Groups of years, shown within ovals are more similar (>65% Bray-Curtis similarity) to ring-mates than to nearest neighbors.Species and corresponding colors are organized by taxonomic groups; mollusks are represented in blue shades, crustaceans in gold and brown, annelids are represented by polychaetes in purples, oligochaetes in greens, and other species are in remaining colors.Wet and dry years are indicated by "W"and "D," respectively.
The upper San Francisco Estuary with four monitoring sites indicated in red Figure 1

Table 1
Species contributing to dissimilarity in assemblage composition among different year categories.Column headers refer to particular combinations of pre-vs.post-Corbula invasion and wet vs. dry year categories.The mean and standard error (SE) of organism abundance per m 2 were computed by species from years in each combined category.Bold values indicate organisms that contributed >2% to dissimilarity in assemblage structure between wet and dry years.Underlined values indicate organisms that contributed >2% to dissimilarity in assemblage composition pre-versus post-Corbula invasion.The species lists are based on similarity percentage analyses presented in Appendices 3 and 4. The SPB site has no pre-invasion data.