Structurally complex habitats, such as seagrasses, provide nurseries that enhance the survival of coastal marine fishes and invertebrates during their early life (Thayer et al. 1984; Bell and Pollard 1989; Gillanders 2006). Investigations of fish communities associated with seagrass beds along the western Atlantic Ocean (Adams 1976; Wyda et al. 2002; Heck and Orth 2006) and other parts of the world (Bell and Pollard 1989; Tolan et al. 1997; Baden and Bostr¨om 2001) document the attributes that seagrasses provide as nursery habitats (Heck et al. 2003). These include refuges from predation, breeding areas, enhanced prey availability, and improved water quality, thereby demonstrating their importance as productive and stabilizing components of the marine environment (Orth et al. 2006).

However, seagrass habitats have been experiencing worldwide declines via escalating threats from anthropogenic influences including direct and indirect effects of chemical pollutants (i.e., nutrient enrichment, contamination) and increasing sedimentation (Ralph et al. 2006; Waycott et al. 2009). Global warming may also alter seagrass species composition by eliminating or displacing species intolerant of warming temperatures or through extreme climatic events (Duarte et al. 2006;Waycott et al. 2009; Diez et al. 2012). These threats endanger not only the seagrasses, but also the associated fish species that rely on these habitats.

Numerous investigations have quantified fish associations and changes in assemblages within seagrass habitats. The most often cited factors affecting fish assemblages include feeding behavior (Grenouillet and Pont 2001; Nagelkerken et al. 2006) and physical gradients (Grenouillet and Pont 2001; Grubbs and Musick 2007). Many investigations were conducted over broad spatial areas but were temporally constrained (Bloomfield and Gillanders 2005; Franca et al. 2009; Pereira et al. 2010; Gray et al. 2011),whereas others have been temporally robust but spatially limited (Fodrie et al. 2010; Sheppard et al. 2011). A study that compared fish communities sampled in 1970 (Livingston 1982, 1985) to fish assemblages in 2006–2007 demonstrated a poleward shift of 13 species indicative of range expansion due to global temperature change (Fodrie et al. 2010). Manipulative experiments in mesocosms have confirmed that species such as Pinfish Lagodon rhomboides and Atlantic Croaker Micropogonias undulatus choose seagrass habitats based on abiotic factors (dissolved oxygen) coupled with biotic (food availability, predation risk) influences (Froeschke and Stunz 2012). These studies document the reduced juvenile fish survival and altered species composition in seagrass habitats that favor warmwater species assemblages due to impacts of anthropogenic stressors and climate warming.

Concern exists in Chesapeake Bay, theworld's second largest estuary, over the decline of seagrass beds since the 1960s, caused principally by light attenuation due to elevated anthropogenic inputs of sediments and nutrients (Orth and Moore 1983; Kemp et al. 2005; Orth et al. 2010). The effects this decline may have had on associated fish fauna, particularly those of commercial or recreational importance, remain poorly documented. Most studies of Chesapeake Bay habitats have focused on single species (Dorval et al. 2005b, 2007; Grubbs and Musick 2007), on a few species (Woodland and Secor 2011), or on lower trophic levels (Kimmel et al. 2006). Although there are valuable studies of commercially important juvenile fish—habitat relationships in Chesapeake Bay, e.g. Atlantic menhaden (Love et al. 2006), few (Orth and Heck 1980; Heck and Thoman 1981; Sobocinski et al. 2013) have examined assemblages associated with seagrass beds. Those that have sampled fish on seagrass beds have posed single-species hypotheses (Dorval et al. 2005b, 2007; Smith et al. 2008) related to growth processes rather than teasing apart the potential factors affecting multispecies juvenilefish assemblages in these habitats or have examined community structure on a limited geographic scale (Orth and Heck 1980; Heck and Thoman 1981). Fishes in ChesapeakeBay use seagrass beds seasonally with the greatest densities of young-of-the-year fish occurring in submerged aquatic vegetation (SAV) from late spring through fall (Orth and Heck 1980; Chesapeake Executive Council 1990). Overall, studies on fish distributions in seagrass habitats throughout the bay are limited and no synoptic investigations exist on fish associations within SAV beds on a broad geographic scale over several years.

From this multiyear study (1997–1999),we provided a broadscale, synoptic evaluation of seagrass-associated fish communities in all major SAV habitats throughout the lower Chesapeake Bay.We examined the effects of physical (salinity, temperature), biological (presence of macroalgae), geographical (zone), and temporal (year) factors on fish abundance within these seagrass beds and tested the null hypothesis of a random fish distribution throughout lower Chesapeake Bay. We also compared the fish community from our collections to historical (Orth and Heck 1980;Weinstein and Brooks 1983) and contemporary (Sobocinski et al. 2013) collections to make inferences about community structure over time.

STUDY SITES

All sites we sampled were characterized by mixed beds of eelgrass Zostera marina and widgeongrass Ruppia maritime (Orth and Moore 1988). Fish species were sampled in the polyhaline—mesohaline lower portion of Chesapeake Bay SAV beds at random locations (Figure 1) nested within three distinct zones (Dorval et al. 2005a, 2005b, 2007; Hannigan et al. 2010). Zone 1 included Tangier and Smith Island in the midportion of the bay (including Bloodsworth Island in 1999); Zone 2 comprised the eastern shore from Crisfield, Maryland, to Cape Charles, Virginia; and Zone 3 encompassed the western shore with its northern boundary at either the Rappahannock River (1997) or Great Wicomico River (1998 and 1999) and southern boundary at Back River. Across these zones there were no differences in seagrass bed density, size, or species composition (Orth et al. 1996, 1997, 1998). These zones are spatially separated by large, deep expanses of the estuary (i.e., river mouths) that likely prevent cross-zone fishmovements (e.g., Dorval et al. 2005b), thereby maintaining the integrity of fish communities on small spatial scales during nonmigratory periods (i.e., summer months).

METHODS

Diurnal, bay-wide fish community sampling was conducted once in 1997 (September) and twice in 1998 and 1999 (August and September). Each synoptic survey took place over 4–5 d during periods of expected high juvenile fish abundance (Orth and Heck 1980). A 4.9-m-wide otter trawl with a 12.7-mm stretch-mesh net and a 6.4-mm stretch-mesh cod end liner was towed at a standardized speed of 1,200 rpm (approximately 4.8 K/h), resulting in a similar area swept at each station. An average water depth of 0.6 m at mean low water was required for trawling effectively during various tidal stages. Trawls were conducted ±2 h of high tide to minimize tidal-related impacts to fish community structure. Fishes were processed onboard immediately after collection, counted as numbers per individual species and returned to the water to minimize the impact on community structure. Salinity and temperature were quantified with a refractometer and stem thermometer, respectively, for most stations. In some cases, if a station was within 500 m of another station (13% of all stations), we assumed that physical characteristics were similar and the salinity and temperature measurements from the nearby station were used. In very few cases (<1% of all stations), temperature and salinity were not directly quantified for a station or a nearby station. For these locations, contemporaneous physical data from the nearest Chesapeake Bay Program monitoring stations ( www.chesapeakebay.net/datawaterquality.aspx) were used.

FIGURE 1.

Map of lower Chesapeake Bay with zones and a typical array of sampling stations (August 1999) indicated.

Seagrass beds were sampled in proportion to their areal coverage determined from the previous year based on annual SAV distribution and abundance mapping efforts (Orth et al. 1996, 1997, 1998). At least 21 sampling locations were randomly assigned within each zone for each year in beds designated as having a cover density of 70–100%, determined from the mapping efforts, which resulted in 545 tows over the course of this study (Table 1). After a tow, we noted the presence of macroalgae in the sample.

In this study, we sequentially sampled the demersal, mobile component of the fish community.We intentionally disregarded more sedentary species (syngnathids, gobies [family Gobiidae], and blennies [family Blenniidae]) because of gear escapement or potential sampling bias for these species. In contrast to other species in the community, their size or close association with the bottom below the seagrass canopy would compromise relative abundance estimates of those individuals. Thus, we are unable to assess the relative abundance estimates of sygnathids, gobies, or blennies in our study relative to those reported in the literature for Chesapeake Bay (Orth and Heck 1980; Weinstein and Brooks 1983), and we have recalculated species abundance after excluding these species for historical comparisons.

TABLE 1.

Number of sites sampled (N), mean temperature and salinity, and percent of sites with macroalgae present for each zone and year sampled.

Statistical analyses.—Redundancy analysis (RDA) was applied to species and environmental data matrices to reveal plausible relationships. Redundancy analysis is a constrained ordination method that models the response (i.e., species matrix) variables as a function of the explanatory (i.e., environmental matrix) variables (ter Braak 1986; Legendre and Legendre 1998). The ordination finds the combination of variables that best explain the variation of the response variables and uses Monte Carlo permutation tests to determine the statistical significance of themodel and each of the explanatory variables. The major advantages to using ordination methods for multivariate data are that transformations are not necessary to fulfill statistical assumptions because statistical significance is assessed with randomization tests and relationships between the response and explanatory data matrices are easily visualized with biplots.

Biplots were constructed with explanatory variables plotted as vectors (continuous) or centroids (discrete), where the vector lengths indicated the relative strength of the relationship with the response data. Response variables are typically plotted as points so that the strength of their relationship with the measured explanatory variables can be visually assessed in the multivariate space by the biplot. Angles between the response variables (plotted as a vector) and explanatory vectors reflect their correlations: correlation is positive when the angle is less than ±90 degrees; correlation is negative when the angle is greater than ±90 degrees.

We fit a model where the species matrix was a function of the environmental parameters (salinity, temperature, presence of macroalgae, year, and location). Salinity and temperature are continuous variables and were plotted as vectors along with the presence of macroalgae.We used an effects model and included year (1997 = 0, 1 or 1998 = 0, 1) and location (Zone 1 = 0, 1 or Zone 2 = 0, 1) as indicator variables. Only two indicators are needed for both the 3 years and three locations to prevent multicollinearity, but all are plotted for clarity. We tested each parameter in the model with 10,000 permutations.

The RDA results were then used as an exploratory analysis prior to generalized additive modeling (GAM). The GAMs are very flexible, and provide an excellent fit when nonlinear relationships and significant noise occur in the predictor variables (Hastie and Tibshirani 1990). Binomial GAMs were developed for each species in the assemblage that showed a relationship with a measured gradient. Local occurrence (presence = 1 and absence = 0) was modeled against environmental variables for all zones. We used a nominal α = 0.05 to assess statistical significance. All statistical tests were conducted with the computer program R (R Development Core Team 2005).

RESULTS

We collected 21,153 fish representing 31 species (Table 2) over the 3 years from 545 otter trawl samples. The overwhelming majority of all individuals collected were Silver Perch (86.1%). Spot made up 5.4% of all individuals collected followed by Weakfish (2.4%), Spotted Seatrout (1.2%), Atlantic Croaker (0.9%), and Atlantic Spadefish (0.9%). The remaining 26 species collectively made up less than 3.0% of all individuals collected (Table 2). We also found that the species composition has shifted from a Spot-dominated community in the late 1970s to early 1980s (Orth and Heck 1980; Weinstein and Brooks 1983) to a community dominated by Silver Perch.

We saw a similar pattern with site occupancy, where the numerically dominant species were found at the greatest proportion of sites (Table 2). Silver Perch were found at >75% of all sites followed by Spot at about 50% of all sites. Weakfish and Spotted Seatrout occurred in about 25% of all sites, while Atlantic Croaker, Southern Kingfish, Summer Flounder, and Atlantic Spadefish occupied approximately 15% of all sites. No other species occurred in more than 10% of sites sampled.

TABLE 2.

Common scientific names of species captured in Chesapeake Bay, frequency of occurrence (%) among all fish captured, and frequency of occurrence (%) among all sites sampled.

There were significant differences between environmental variables among zones (Table 1; F_6,1078 = 32.58, P < 0.0001). Temperature (F_3,540 = 14.34, P < 0.0001), salinity (F_3,540 = 78.20, P < 0.0001), and the proportion of sites with macroalgae present (F_3,540 = 3.35, P = 0. 0187) showed differences among zones. Year was a significant covariate only for salinity (F₁ =76.09, P <0.0001). Salinitywas different among all zones and was higher during 1999 than other years. Temperature was highest in zone 3 and similar between zones 1 and 2. There were more sites in zone 3 containing macroalgae than at either zone 1 or 2, where proportions were similar (Table 1).

TABLE 3.

Redundancy analysis eigenvalues, cumulative percent of variance explained, and Monte Carlo permutation tests (10,000 permutations) for all axes.

The first four axes of the RDA ordination were significant (Table 3) and explained91%of the cumulative variance. The first axis explained 53% of the total variance. All measured explanatory variables explained a significant amount of the variance (Table 4). Salinity appeared to be the most important variable structuring the fish community (Figure 2). Temperature and the indicator variables zone 1 and 1997 were also important. Zone 1 was negatively correlated with salinity and four species of fish (Weakfish, Atlantic Croaker, Spotted Seatrout, and Southern Kingfish) were most often encountered in this sampling area. Conversely, Spot occurred most often at higher salinity sites in zone 3 and Silver Perch were associated with moderate to high salinity sites in zone 2.

Because all effects in the RDA model were significant, we built GAMs for the 10 species that showed a relationship with measured gradients (Table 5). The deviance explained for these models ranged from 3.2% for Spotted Seatrout to 35.1% for Pigfish. Both temperature and salinity were significant for seven or six species, respectively. Most species show a negative or no response to low temperatures while showing a variable response at moderate to high temperatures (Figure 3). For salinity there was a mixed response at both high and low salinities (Figure 4). Using 1999 as the baseline for comparisons, 1998 differed the most as five species were either more (three) or less (two) likely to be present at sampling locations. Conversely, during 1997 only Spot were significantly more likely to occur at sampling locations. Similarly, using zone 3 as a baseline, five species were either more or less abundant in at least one of the other two zones. Weakfish were much more likely to be present in both zones 1 and 2 than zone 3, while Pigfish were much less likely to be present in zones 1 and 2 than zone 3. The other species showed a mixed response. The presence of macroalgae only negatively impacted the presence of Spot.

TABLE 4.

Monte Carlo permutation tests (10,000) for each term in the RDA model. The terms 1999 and zone 3 were not included in the model because of multicollinearity.

DISCUSSION

Juvenile finfish populate nursery grounds in Chesapeake Bay in response to habitat complexity (Orth and Heck 1980; Weinstein and Brooks 1983), quality (presence of macroalgae; Sogard and Able 1991), and gradients of temperature and salinity. We explicitly tested the null hypothesis that the fish community would not respond to abiotic gradients but found that responses were highly variable between species where site occupancy was influenced by the temperature and salinity regime. For example, we showed that Spot responded positively to both warmer temperatures and higher salinities. Under a climate change scenario that included increased precipitation, Spot would not be favored. In contrast to our findings, previous research in Chesapeake Bay and the Hudson River estuary indicated that the juvenile finfish community largely responded to abiotic conditions over the previous year and current conditions had no effect on community dynamics (Hurst et al. 2004; Wingate and Secor 2008). Our research demonstrated the advantage of multiple-year broad-scale synoptic sampling of these nursery areas by quantifying potential abiotic drivers of community structure in the current year.

FIGURE 2.

Redundancy analysis biplot examining juvenile finfish occurrence in Chesapeake Bay in relation to measured environmental parameters. Continuous response variables are represented by vectors and dummy indicator variables are represented by a plus sign for year (1997–1998) and a multiplication sign for location (Zone 1–3). See Table 2 for species codes.

TABLE 5.

Generalized additive modeling results for the effects of environmental covariates on species presence—absence. Percent deviance is the percentage of the deviance explained for the model, temperature and salinity are the estimated degrees of freedom of the smoothing parameter, and year (1997, 1998), zone (1, 2), and algae (presence of macroalgae) are the estimated parameter effects for the remaining terms in the model. Values in bold italics indicate significance at α = 0.05 and blank cells indicate that there were no individuals collected during 1997 and therefore only 1 year is needed for an effects parameterization.

The physical and chemical structure of Chesapeake Bay is well known (Austin 2004; Dorval et al. 2005a; Hannigan et al. 2010); there is a salinity gradient along the main-stem portion of the estuary and differences in salinity between areas on the eastern and western shores. The salinity structure is highly dependent on precipitation and discharge from the many tributaries supplying freshwater. Wingate and Secor (2008) have shown that both winter temperature and flow (a proxy for salinity) are drivers structuring the fish communities in the upper portion of Chesapeake Bay. Interestingly, Austin (2004) found that salinity was lagged by about 90 d in response to freshwater influx. This is a similar response to that noted for the fish community in Chesapeake Bay (Wingate and Secor 2008). Likewise, we have demonstrated the importance of salinity for structuring the seagrass-associated juvenile fish communities. Perturbations to the salinity regime in Chesapeake Bay have the potential to alter the value of these nursery habitats to overall stock dynamics. One of the looming drivers of the physical and biological structure of the ecosystem is that from climate change (Pyke et al. 2008). The predictions for the mid-Atlantic region are for increased variability in precipitation, which, as we have demonstrated, will lead to variable responses from the juvenile fish community. The approach we took to use GAM sets up a framework that enabled us to make inferences about how this fish community would respond to potential threats from climate change. We demonstrated that numerous fishes in the community would respond negatively to decreased salinity.

Presence of macroalgae was a significant variable in the analysis for Spot and the overall fish community. Macroalgae generally occur as drift accumulations throughout Chesapeake Bay seagrass beds because seagrasses reduce current speed and buffer wave action. At very high biomass levels, macroalgae can smother and eliminate seagrass (Hauxwell et al. 2001) and also promote hypoxia that negatively influences fish and invertebrate populations (Baden et al. 1990; Oesterling and Pihl 2001; Deegan et al. 2002; Fox et al. 2009). Alternatively, at low to moderate levels of abundance, macroalgae can increase habitat complexity of seagrass habitats and provide additional food resources for resident invertebrates, thereby providing added forage for fish populations (Martin-Smith 1993; Norkko et al. 2000; Epifanio et al. 2003; Powers et al. 2007). The significant negative association of macroalgae with the presence of Spot and the overall fish community suggested that macroalgae was, and will continue to be, a concern in Chesapeake Bay and potentially other estuaries that are or are becoming more eutrophic, which will influence macroalgal abundances (McGlathery et al. 2007).

Two sciaenids, Silver Perch and Spot, were the numerically dominant species in fish assemblages associated with Chesapeake Bay seagrass beds in our study. In this study, Silver Perch were far more abundant than Spot and present at most sampling locations. This pattern contrasts sharply with a study by Orth and Heck (1980), who sampled at similar times using similar gear and reported that Spot were also ubiquitous but they dominated the relative abundance of species encountered in western shore locations (Mobjack Bay and the York River) of Chesapeake Bay from 1976 to 1977 comprising over 63% of individuals collected. Silver Perch in their study (Orth and Heck 1980) made up just over 5%of the fish assemblage. A similar abundance relationship was found byWeinstein and Brooks (1983) within a seagrass bed at the mouth of Hungars Creek on the eastern shore of the bay. They reported that Spot comprised >80% of individuals collected, whereas Silver Perch accounted for between only 1% and 2% of the fishes sampled in seagrass. More recent sampling in seagrass beds in the western shore of lower Chesapeake Bay between 2009 and 2011 demonstrated continued dominance of Silver Perch (Sobocinski et al. 2013), suggesting a long-term, dramatic reversal in relative abundance of the two species.

FIGURE 3.

Generalized additive modeling results for the effect of temperature on the presence of selected species that showed a relationship with a measured environmental gradient. The statistical significance of the smoothing function is indicated in Table 5. The dashed lines around the smoothed fit are 95% confidence intervals, and data availability is indicated by tic marks above the x-axis. See Table 2 for species codes.

Other species frequently encountered in our study included Weakfish, Spotted Seatrout, Northern Kingfish, Summer Flounder, Atlantic Croaker, Atlantic Spadefish, and Black Sea Bass. Several of these (Northern Kingfish, Atlantic Spadefish, and Weakfish) were not collected in earlier investigations (Orth and Heck 1980;Weinstein and Brooks 1983) and there may be multiple reasons why these fish were not observed. First, the lack ofWeakfish in the historic surveys likely represents sampling in seagrass habitats that are not preferred.We found thatWeakfish were found at very few sites on the western shore (zone 3) or eastern shore (zone 2), which corresponds to areas sampled in historic surveys (Orth and Heck 1980; Weinstein and Brooks 1983), but site occupancy was generally high in zone 1. Second, Atlantic Spadefish were not collected in historic surveys but were found in at least 14% of all sites we sampled during 1998–1999. However, Atlantic Spadefish were not captured at any site we sampled during 1997. Therefore, it is possible that Atlantic Spadefish were not present during the years sampling took place in the historical surveys. Finally, the current presence of Northern Kingfish in seagrass beds may also represents an apparent change in species composition over time. Northern Kingfish were present in all geographic zones and occupied between 7% and 33% of all sites sampled. Any of these three reasons (unfavorable habitat, poor recruitment year, and change in distribution) could be applied to any of these three species, and it underscores the importance of our multiyear synopticscale approach to addressing hypotheses concerning community structure. However, these three species also were abundant in the 2009–2011 contemporary study (Sobocinski et al. 2013) again offering evidence of a species reversal in these seagrass fish communities.

FIGURE 4.

Generalized additive modeling results for the effect of salinity on the presence of selected species that showed a relationship with a measured environmental gradient. The statistical significance of the smoothing function is indicated in Table 5. The dashed lines around the smoothed fit are 95% confidence intervals, and data availability is indicated by tic marks above the x-axis. See Table 2 for species codes.

A valuable use for studies like ours when similar historical studies exist is for comparisons of taxa to assess community change in response to climate change (Fodrie et al. 2010). However, the major drawback to these comparisons is that they will involve species that are rare and analyzing data with many zeros presents difficulties (Pennington 1983). Our study location, Chesapeake Bay, is ideally positioned just north of Cape Hatteras, North Carolina. Cape Hatteras has long been recognized as a faunal break (Grothues and Cowen 1999) and stock boundary (Bowen and Avise 1990; Jones and Quattro 1999). Therefore, the impacts of climate warming and range extension of southerly species, as well as range contraction of northerly species, should have the highest likelihood of detection here. As pointed out by Fodrie et al. (2010), an increased abundance of southerly species can be used as evidence for the effects of climate change. Certainly Atlantic Spadefish follow this pattern because Chesapeake Bay is at the northern end of their range and they were present at a moderately high proportion of sites and made up >1% of the total numbers of individuals collected in this study. In the Sobocinski et al. study (2013), Atlantic Spadefish increased further to 10.4% of the collection. Similarly, three other species (Florida Pompano, Gray Snapper, and Mojarra) were also found in our study, although at low abundances and among few sites. Although none of these species were found in historic collections (Orth and Heck 1980; Weinstein and Brooks, 1983), they are a part of the Chesapeake Bay fauna. To be able to conclusively take this as evidence of climate change, we must establish that these juveniles are not vagrants (Sinclair 1988; Sinclair and Iles 1989) and have either established reproducing populations or are contributing to the spawning population. Currently, we do not have information concerning the overwinter survival and contribution of southerly spawned juveniles transported to Chesapeake Bay, but we cannot discount the evidence that seems to indicate that their presence may be increasing. It would likely require a technique such as otolith chemistry (Schaffler et al. 2009) to definitively answer the question of whether these juveniles are contributing to the adult population.