Linking the Abundance of Estuarine Fish and Crustaceans in Nearshore Waters to Shoreline Hardening and Land Cover

Human alteration of land cover (e.g., urban and agricultural land use) and shoreline hardening (e.g., bulkheading and rip rap revetment) are intensifying due to increasing human populations and sea level rise. Fishes and crustaceans that are ecologically and economically valuable to coastal systems may be affected by these changes, but direct links between these stressors and faunal populations have been elusive at large spatial scales. We examined nearshore abundance patterns of 15 common taxa across gradients of urban and agricultural land cover as well as wetland and hardened shoreline in tributary subestuaries of the Chesapeake Bay and Delaware Coastal Bays. We used a comprehensive landscape-scale study design that included 587 sites in 39 subestuaries. Our analyses indicate shoreline hardening has predominantly negative effects on estuarine fauna in water directly adjacent to the hardened shoreline and at the larger system-scale as cumulative hardened shoreline increased in the subestuary. In contrast, abundances of 12 of 15 species increased with the proportion of shoreline comprised of wetlands. Abundances of several species were also significantly related to watershed cropland cover, submerged aquatic vegetation, and total nitrogen, suggesting land-use-mediated effects on prey and refuge habitat. Specifically, abundances of four bottom-oriented species were negatively related to cropland cover, which is correlated with elevated nitrogen and reduced submerged and wetland vegetation in the receiving subestuary. These empirical relationships raise important considerations for conservation and management strategies in coastal environments.


Authors Authors
This article is available at W&M ScholarWorks: https://scholarworks.wm.edu/vimsarticles/763 Introduction Coastal systems provide important ecosystem services with considerable economic and ecological value (Lotze et al. 2006; Barbier et al. 2011). Thus, it is critical that we understand how accelerating land-based stressors will affect living resources in coastal ecosystems. Estuaries and other coastal areas are rapidly changing in response to intensifying stressors linked to human activities at both local and global scales (Lotze et al. 2006;Halpern et al. 2008;Barbier et al. 2011;Allan et al. 2013). Human population centers are already located predominantly in coastal areas, and population density near coasts continues to rapidly increase (Crossett et al. 2004). Sea levels are projected to rise substantially by mid-century (IPCC 2007), which will amplify demand for shoreline hardening as humans seek to protect cities, residences, and agricultural lands from coastal flooding (Rahmstorf 2007;Arkema et al. 2013). The ability to conserve coastal ecosystems is further complicated by uncertainties associated with predicting how humans and ecosystems will respond to these and other elements of global change (Runting et al. 2013).
Human population growth will necessitate increased urban and agricultural land use and without careful management can be expected to accelerate land-based pressures such as eutrophication and shoreline hardening on coastal systems. Urban and agricultural lands contribute to increased loads of nutrients, contaminants, and sediments; rapid changes in surface runoff and river discharge; and shoreline habitat alteration (Foley et al. 2005). Nutrient loading resulting from landbased activities can reduce water clarity, dissolved oxygen, and submerged aquatic vegetation (SAV), increase harmful algal blooms, and alter food webs (Nixon 1995, Cloern 2001Kemp et al. 2005). Elevated concentrations of metal and organic contaminants can accumulate in tissues of aquatic organisms, affecting organism populations and human health (King et al. 2004). Impervious surfaces exacerbate runoff of nutrients, contaminants, and fine sediments, and contribute to rapid changes in runoff and water discharge (Wheeler et al. 2005;Uphoff et al., 2011). Thus, land use change driven by human population growth has the potential to affect fish and crustaceans (Caddy 1993;de Leiva Moreno et al., 2000;Breitburg et al. 2009;). Unfortunately, direct links between land cover and fish and crustacean population changes have been difficult to identify because of challenges in interpreting fisheries landings data (de Mutsert et al. 2008) and in scaling up from experiments and local observations to larger spatial extents (Boynton et al. 2001;Breitburg et al. 2009).
Shoreline habitat alteration, including replacement of natural shorelines with hardened structures, is also growing at an alarming rate in coastal areas (Brody et al. 2008;Gittman et al. 2015) and will accelerate with global sea level rise (Rahmstorf 2007). Natural habitats including shoreline vegetation and intact reefs offer some protection against sea level rise, but evaluations of these natural defenses lag behind those for shoreline hardening (Jones et al. 2012;Arkema et al. 2013). Critical habitats have already been subjected to substantial losses, as an estimated 50-65% of tidal wetlands have been lost or degraded worldwide (Lotze et al. 2006;Gedan et al. 2009; Barbier et al. 2011). Shoreline hardening degrades environments by reducing water and sediment quality, altering communities of benthic infauna, and contributing to loss of emergent wetland and SAV (Peterson & Lowe 2009). However, the effect of shoreline hardening on assemblages of fishes and crustaceans is poorly understood (Seitz et al. 2006;Bilkovic & Roggero 2008;Strayer et al. 2012). Localized negative effects have been documented for early life stages (Rozas et al. 2007) and in areas where hardening alters the subtidal zone (Toft et al. 2007), which can affect the abundances of fish species that select for specific substrates (Munsch et al. 2015a). Cumulative, system-scale impacts of shoreline alteration on aquatic fauna may be substantial (Peterson & Lowe 2009;Dethier et al. 2016). Unfortunately, systemscale effects are not easily identified from local-scale studies due to difficulties in scaling up or in designing a study with sample sizes adequate for capturing the distribution of the response variable (Breitburg & Riedel, 2005;Schlacher et al. 2008). What is really needed is a meta-analysis approach that captures the variation that exists in land-use and shoreline hardening with enough sampling effort to detect their potential effects on aquatic systems at larger spatial scales.
Concurrent examination of watershed land use and shoreline habitat alteration at landscape scales is especially important (Cloern 2001;Bilkovic & Roggero 2008) because of the synergistic effects of multiple stressors on coastal environments (Breitburg & Riedel 2005;Lotze et al. 2006;Halpern et al. 2008). Biota from shallow nearshore zones (highly productive interfaces between land and water [Beck et al. 2001]) may show the strongest response to land-based anthropogenic stressors, and if so they could thus serve as indicators of their ecological effects (Diaz & Rosenberg 2008;Strayer et al. 2012). Shallow water nearshore zones directly interact with terrestrial ecosystems, are structurally altered by shoreline hardening, and are the first areas to receive effluent and runoff from upland areas. Comprehensive landscape-scale analyses addressing the effects of anthropogenic stressors on mobile fish and crustaceans are needed to identify causes of impacts, and to enable estimation and public education of benefits of conservation and restoration initiatives.
Here, we present a cross-system comparison examining the relationships between mobile fish and crustacean abundances, land cover (e.g., developed [i.e., areas with constructed materials and impervious surfaces, such as urban, suburban, industrial or commercial areas] and agricultural land vs. forest and wetland), and localand system-scale shoreline habitat alteration. We focus on subestuaries (i.e., tributary estuaries) of Chesapeake Bay and Indian River Bay (one of the Delaware Coastal Bays). Like other coastal areas, these systems provide an enormous economic benefit. Chesapeake Bay, for example, supports hundreds of thousands of jobs and provides an annual economic benefit of 60 billion USD (Chesapeake Bay Blue Ribbon Finance Panel 2004). Our dataset (Fig. 1) includes over 600,000 individuals sampled from 587 sites within 39 subestuarine systems, providing a powerful basis for inferring the relative impacts of multiple stressors on the abundance of estuarine fishes and crustaceans (i.e., crabs and shrimp). Our approach (1) identified key local-scale differences in communities of fish and crustaceans adjacent to natural vs. hardened shorelines and (2) demonstrated that systemscale faunal abundance is often related to cumulative shoreline condition (defined in this paper as the percentage of system's shoreline comprised of hardened and natural shoreline types) and/or watershed cropland cover. This study establishes that mobile fish and crustacean abundances are tightly linked with multiple stressors related to human activities, including shoreline hardening and watershed land use, and has implications for the study, conservation, and management of coastal living resources under global change.

Methods
We examined patterns in the abundances of nearshore estuarine fauna from 587 sites and 39 subestuaries within Chesapeake Bay and the Delaware Coastal Bays with surface areas from 1.8 to 100.8 km 2 , and salinities of 6.0 to 20.0 PSU (Fig. 1, Online Resource A [ Table A1]). Subestuaries were deliberately chosen to span broad gradients of crop and developed land cover across watersheds National Land Cover Database, Fry et al. 2011. This design enabled a cross-system comparison approach in which subestuaries/coastal bays were the units of replication for evaluating relationships of faunal abundances with watershed land cover and cumulative shoreline condition. Fish and crustacean abundance data were from a composite of six contributed datasets using either beach seines or fyke nets to examine nearshore fish and crustacean communities (see Online Resource A [ Table A2] for a description of sampling gears and procedures from each of the six contributed datasets). An objective of each of these studies was to evaluate linkages between land cover patterns and mobile estuarine fauna in nearshore waters. Each contributing dataset used a unique method for sampling fauna, and data were collected over multiple years (2002-2003 and 2008-2012). In addition, the efficiency and catchability of seine nets were likely to differ among the four shoreline types assessed. Therefore, the following are the three major challenges that needed to be addressed to analyze this composite dataset: (1) probable differences in beach seine efficiency at different shoreline types, (2) potential differences in sampling efficiency across gear types used, and (3) likelihood of interannual variability in fish and crustacean abundance.
Our methods accounted for differences in sampling efficiency among gear types, sampling years, and shoreline types, and their interactions, by using Leslie depletions and z-scoring (Online Resource A). Briefly, abundances were adjusted for shoreline-specific catchability bias, as determined by 4× Leslie depletions conducted at 36 sites (9 each of beach, bulkhead, riprap, and wetland, Online Resource A [ Table A3]), and then log-transformed so there were no glaring violations of the assumption of normality. Abundance data were then zscored within each unique combination of gear type and year of capture (total of 14 such combinations) to put abundance distributions from different collection methods on the same scale and to remove variation attributable to sampling methodology and inter-annual variability. Neither gear type nor year of capture was a significant predictor of z-scored abundance or zscored fish richness (two-factor ANOVAs; all P > 0.80 for both factors). Z-scoring produced the following index of species abundance: where Abundance average is the average catch of a particular gear in a single year and σ is the standard deviation of abundance values from that gear/year. This procedure converted abundance values to standard deviation units with respect to the average catch of a gear type in a particular year. Z-scores (or Bresponse ratios^) provide a comparison of relative rather than absolute abundance and are commonly used in meta-analyses for comparing data collected using different protocols (e.g., Felker-Quinn et al. 2013;Shade et al. 2013). See Online Resource A for additional details on treatment of faunal abundance data. Fig. 1 Map of sample locations and land use. All sites were used for system-scale analyses. Purple diamonds represent sites that were also used in assessing shoreline type patterns within 0-3 and 0-16 m of shore.
Blue squares represent sites that were also used in assessing shoreline type patterns farther from shore (10 to 50 m distance). Black circles represent sites that were only used for system-scale analyses. Samples were collected over multiple years, including 2002-2003 and 2008-2012, by multiple field survey teams from throughout the Chesapeake Bay region. Land use National Land Cover Database, Fry et al. 2011 is summarized as developed land (gray), cropland (orange), forest (dark green), and nearshore wetland (light green). Pasture land (light yellow) is also included for reference, but not included in analyses. See Online Resource A for attributes of each subestuary and details on how land use categories were used in the analysis A total of 15 taxa were analyzed across the different datasets (species names provided in Tables 1 and 2), including 13 species (11 finfish and 2 crustaceans) plus the family Centrarchidae (sunfishes; usually Lepomis gibbosus [pumpkinseed] or L. macrochirus [bluegill]) a n d t h e g e n u s M e n i d i a ( A t l a n t i c s i l v e r s i d e [M. menidia] and inland silverside [M. beryllina]), which are lumped due to difficulty identifying small individuals to species while field-processing large numbers of live fish. Taxa included in analyses were chosen based on prevalence in the study area, such that sample sizes were adequate to examine species abundance patterns with shoreline type and across broad subestuaryscale gradients. We constrained our analysis of abundance patterns to common taxa that occurred in at least 65% of subestuaries and at least 20% of sites (these cutoffs were natural breakpoints in the data), with the exception of fish species richness, which included all taxa. This was done to ensure that the species included in the analysis spanned land use, shoreline hardening, and latitudinal gradients present in the study design. Taxa were divided into three functional groups for the purposes of comparison (Table 1). Analysis was limited to samples collected between June 26 and October 10 to capture relatively stable summertime temperature and fish distribution patterns and to avoid transitional spring and autumn periods. Most of the six contributing datasets did not collect specimens outside of this date range, precluding a formal assessment of seasonal effects. This study includes two main analyses evaluating (1) local-scale patterns (i.e., those in waters directly adjacent to shore) and (2) system-scale patterns (i.e., those within a whole subestuary) in fish and crustacean abundance, as well as several a posteriori analyses to elucidate mechanisms.

Local-Scale Effects of Shoreline Type
We evaluated the local-scale effects of four common shoreline types on fish and crustacean abundances for 12 species. Two natural shoreline types (native wetland and beach) and two hardened shoreline types (bulkhead and riprap revetment) were considered. Local-scale effects of shoreline type were analyzed separately from system-scale effects to capitalize on the strengths of two of the contributing datasets that were explicitly designed to test for effects of shoreline hardening at sampled sites. Although shoreline type was recorded at all sample sites, the other four datasets included in our metaanalysis were not designed to test the effects of shoreline type and had varying numbers of sites from each shoreline type in each subestuary. By contrast, the two datasets (total of 16 subestuaries and 260 samples using beach seines deployed parallel to shore) had a balanced complete block design to test for shoreline type effects, with equal representation of each shoreline type in each subestuary, and with sites in proximity to one another within each subestuary to limit extraneous effects of site location (see Online Resource A for details). These two datasets also explicitly sampled in two different zones of water that may vary in their use by mobile fish and crustaceans. Fifteen-meter seine nets that sampled the area of water extending from the shoreline to 3 m out from shore (hereafter Bwithin 0-3 m^; n = 36 for each habitat type) were used to describe the very narrow band of aquatic habitat at the land-water interface, while 60-m seine nets that sampled from the shoreline to an average distance of 16 m out from shore (hereafter Bwithin 0-16 m^; n = 29 for each habitat type) were used to describe a broader area of nearshore habitat (see Online Resource A for details). Three taxa included in the system-scale analyses (Atlantic croaker [Micropogonias undulatus], Atlantic menhaden (Brevoortia tyrannus), and silver perch [Bairdiella chrysoura], see next section) occurred at Table 1 Generalized functional groupings of studied species, based on descriptions from Murdy et al. (1997) (Masterson 2008) b (Jordan 2009) <15% of sites in the local-effects dataset and were excluded; each of these species had low abundance even where they were present in this portion of the dataset. Single-factor ANOVAs blocked by subestuary of capture were used to assess the effects of shoreline type, a categorical independent variable, on species abundances, fish species richness, and local-scale water depth. Tukey HSD procedures were then used to identify significantly different shoreline type pairs. This analysis method was chosen because the significance testing best identifies species-specific shoreline type patterns, which are likely of interest to readers. However, we also evaluated the local-scale effect of shoreline type on faunal community structure as a whole using redundancy analysis (Borcard et al. 1992), a multivariate approach useful for understanding patterns in community composition that has been frequently applied to fish communities (e.g., Angermeier and Winston 1999;Sharma et al. 2011;Kornis et al. 2013). We used redundancy analysis to determine the total amount of variation in faunal community structure explained by local-scale shoreline type and to corroborate the results of the ANOVA approach (See Online Resource C for more detail).
To augment our analysis of local-scale response of species to shoreline type, we also analyzed a third dataset included in the system-scale analysis (26 subestuaries and 222 samples, fyke nets set parallel to shore and soaked for 24 h [Online Resources A, Table A2]) that sampled faunal abundance at 10 to 50 m distance from shore, where water depths were approximately equal across a gradient of wetland shoreline occupation. This dataset was not designed to test for local effects of shoreline hardening per se, but instead examined the amount of shoreline habitat occupied by emergent wetlands, assessed visually at each site, on a continuous percentage scale. We used this dataset to evaluate species abundance patterns with the amount of emergent wetland (linear regression) in the 10 to 50-m offshore stratum where water depth did not co-vary with shoreline type, thus aiding in the interpretation of shoreline type effects with respect to water depth. The samples from the 10 to 50-m offshore stratum also facilitated ΔAIC scores are in colored boxes and are the difference in AIC between the best model and the model excluding each predictor, thereby approximating each predictor's relative influence (Cohen et al. 2012). ΔAIC scores are not multiplicative (e.g., a score of '6' is not three times greater than a score of B2^). As a rule of thumb, a ΔAIC <2 suggests removing the predictor still produces a model with substantial support; ΔAIC values between 3 and 7 indicate that a model without the predictor has considerably less support; and ΔAIC >10 indicates that dropping the predictor produces a model that is very unlikely (Burnham and Anderson 2002). Negative relationships are red; positive relationships are blue. Predictors included in alternate models (which are within 2 AIC units of the best model and are thus considered to be nearly as likely to represent truth as the best model) are stippled red or blue and do not contain ΔAIC scores. Tot Var and Sub Var are the % of total and among-subestuary variance explained by the best model, respectively. δ Centrarchidae may be planktivores or benthivore-piscivores depending on life stage and/or available food resources. Thus, they are placed in both functional groups. *Cropland explained 27.4% of among-subestuary variation in Atlantic menhaden abundance with a high-influence subestuary removed. # % hardened shoreline noted for B. chrysoura and L. xanthurus because of intriguing univariate patterns, but was not included in alternate models. Crop = % watershed crop land; Dev = % watershed developed land; NS-For = % nearshore forest; NS-Wet = % nearshore wetlands; Hard Shore = % hardened shoreline; Beach Shore = % beach shoreline beach. Some candidate predictors (i.e., subestuary depth, subsetuary area, salinity, distance from the ocean, and local-scale shoreline type) included in the models were omitted from this table for brevity and clarity. See Online Resource B (Tables B4 and B5) for complete details on model output comparison of local-and system-scale effects of wetland habitat on select species.

System-Scale Predictors of Fish and Crustacean Abundance
All six contributing datasets (39 subestuaries, 587 sites) were used to evaluate system-scale (i.e., subestuary-scale) predictors of fish species richness and the abundance of 15 commonly occurring taxa in nearshore waters. Multilevel mixed-effect regression models, whose hierarchical structure accounted for site aggregation, were used to examine patterns in each response variable (Proc Mixed, SAS Version 9.2). Subestuary of capture was included as a random effect to account for the greater potential for similarity among sites from the same subestuary. Between-subestuary fixed effects included cumulative shoreline condition (which is comprised of % hardened shoreline, % beach shoreline, and the % of land within 100 m of shore comprised of either wetlands [hereafter, B% nearshore wetland^] or forest [hereafter, B% nearshore forest^]); watershed land cover (developed land and crop land); and average subestuary water depth, surface area, and surface salinity. Distance between the ocean and each subestuary's mouth (hereafter Bdistance from the ocean^) was included in place of salinity for ocean-spawning species that begin life in or near the ocean shelf and return there to spawn because abundances of such species may decrease with increasing distance from the ocean. The within-subestuary factor of local-scale shoreline type (beach, wetland, bulkhead, or riprap) was also included as a fixed effect due to significant patterns identified by the local-scale analysis. Wetlands and forested land within 100 m of shore, measured at the subestuary level, were included as predictors rather than whole-watershed wetland/forested land because natural land covers (e.g., forest and wetlands) close to the shoreline can buffer against runoff of nutrients and other pollutants (Vernberg 1993) and provide critical habitat (woody debris and intertidal wetland surfaces) for aquatic mobile estuarine fauna (Everett and Ruiz 1993;Bilkovic & Roggero 2008). Models that included this full set of candidate predictors satisfied non-collinearity requirements laid out in Weller et al. (2011), and correlation among all potential predictors was limited, with the absolute value of all Pearson's r < 0.65. Please see Online Resource A for details on how system-scale environmental characteristics (e.g., % watershed land use, % hardened shoreline) were determined. Akaike's information criterion (AIC) was used to compare the relative quality of each model. The best model of each response (i.e., species abundances and fish species richness) was determined separately by removing single predictors from the full model until the model with the lowest adjusted AIC score was identified. Models within 2.0 AIC units of the best model (hereafter, Balternate models^) were considered to be nearly as likely as the best model to accurately describe the variability in the response variable (Burnham and Anderson 2002) and thus were also reported. We used a Bdrop-oneâ nalysis in which the difference in AIC between the most likely model and the model minus each individual predictor (hereafter BΔAIC^) is used to evaluate each predictor's relative importance (Cohen et al. 2012). ΔAIC values used in this way reflect the reduction in model likelihood absent a given predictor and, thus, helps determine how much less explanatory power a model has when a given predictor is removed. As a rule of thumb, a ΔAIC <2 suggests removing the predictor still produces a model with substantial support; ΔAIC values between 3 and 7 indicate that a model without the predictor has considerably less support; and ΔAIC >10 indicates that a model without the predictor is very unlikely (Burnham and Anderson 2002). In addition, we report the total variation explained by the fixed effects, as well as the % of variance among subestuaries explained by each model. The following procedures are analogous to calculating marginal R 2 values sensu Nakagawa & Schielzeth (2012). The total variation explained by each model was calculated as: where V best is the residual variance in the best model and V int is the residual variance in an intercept only model (no fixed or random effects). The total variance that occurs between subestuaries (V sub ) was determined as the amount of variance attributable to the random effect of subestuary of capture in a model that included no other effects. The % of among-subestuary variation explained by fixed effects in the best model (B% Sub. Var.^) was calculated as: where V sub-only is the variance explained by subestuary of capture in a model with no fixed effects and V sub-best is the variance explained by subestuary of capture in the best model. This is a conservative evaluation of the importance of the fixed effects because it accounts for variation explained by subestuary of capture prior to determining the amount of variation attributable to fixed effects.
For each species, faunal abundances were modeled on two subsets of the overall dataset to isolate relationships between faunal abundances and either cropland (a Bcropland subsett hat excluded watersheds with >25% developed land cover) or developed land (a Bdeveloped land^subset that excluded watersheds with >18% cropland cover). This was necessitated by a wedge-shaped relationship between cropland and developed land (Fig. 2) where low proportions of crop land were associated with both low and high values of developed land and vice versa. Exclusion cutoffs were determined by visual inspection of apparent thresholds in the relationship between cropland and developed land (Fig. 2). Subsetting is not ideal, and we initially used a principal components analysis (PCA) to avoid subsetting by loading land use onto orthogonal axes. However, the resultant PCA loadings, when used as predictors of species abundances, severely restricted clear interpretation of land use patterns; thus, we chose the subset approach to facilitate interpretation. In the absence of subsetting, low values of one human land use type would be associated with either high amounts of the other human land use type or high amounts of wild land use, which is problematic for detecting patterns of species abundance along human land use gradients. Simply lumping cropland and developed land into one category would also be problematic because cropland vs. developed land effects may be related to different mechanisms. Thus, subsetting prevented possible masking of effects by ensuring that low values of the human land use type (either cropland or developed land) corresponded only with high values of wild land use. In effect, watersheds with low percentages of human land use (i.e., watersheds with both <18% cropland and <25% developed land) became part of a de facto control group of relatively wild watersheds present in both subsets. On average, subestuaries in the Bcontrol group^had 68.0 ± 3.6% (±SE) forested land cover vs.
42.0 ± 3.2% (±SE) in watersheds with >18% cropland and 43.6 ± 4.4% (±SE) in watersheds with >25% developed land. Thus, negative cropland or developed land use effects could also be interpreted as positive effects of wild or forested lands. This is useful in our study because % watershed forest land could not be explicitly considered in our models due to negative collinearity with cropland (R 2 = 0.69, p < 0.0001) and developed land (R 2 = 0.74, p < 0.0001).

A Posteriori Analyses
Several a posteriori analyses were conducted after the main two analyses to help elucidate mechanisms behind identified relationships and to provide supporting information for the interpretation of our findings.
To evaluate competing mechanisms for observed system-scale patterns with watershed cropland and with cumulative shoreline condition, we regressed watershed cropland, watershed developed land, and % hardened shoreline against mechanistic responses including nutrient enrichment (measured as total nitrogen), dissolved oxygen, submerged aquatic vegetation (hereafter BSAV^), and nearshore wetlands. This was done to establish how potential mechanisms correlated with subsetuary-scale land use and shoreline hardening predictors. We used nitrogen to represent nutrient enrichment because it is considered the limiting nutrient and primary driver of eutrophication in coastal systems during summer (Howarth & Marino 2006) and because nitrogen reduction is a common goal of coastal management plans (USEPA 2010).
We further evaluated possible mechanisms behind subestuary-scale patterns between species abundances and watershed cropland by examining the performance of models in which total nitrogen (TN), SAV, both, or both plus their interaction were substituted for watershed cropland in multivariable mixed effect models (described above). If substituting a mechanistic predictor for % watershed cropland improved model performance, it would lend support for that predictor as a mechanism behind relationship between cropland and species abundance. We also determined the % of among-subestuary variation explained by univariate models containing only % watershed cropland, TN, SAV, % nearshore wetlands, or salinity/distance from the ocean. As with the % watershed cropland analysis, only subestuaries with watersheds comprised of <25% developed land were considered (n = 19).
To test whether offshore species abundances were also related to % watershed cropland, we evaluated catch of species in otter trawl samples from offshore mid-channel bottom habitat (n = 9 subestuaries, Online Resource A). This analysis was limited to species that were also Fig. 2 Wedge-shaped relationship between the proportion of crop and developed land within the watershed. Each data point represents one subestuary. Subestuaries were split into two datasets emphasizing gradients of cropland and developed land in order to avoid potential confounding effects of each anthropogenic land use type. Apparent thresholds were used to determine subset composition. To examine species relationships with cropland, subestuaries with >25% developed land (points above the horizontal dashed line) were excluded. Similarly, to examine species relationships with developed land, subestuaries with >18% crop land (points to the right of the vertical dashed line) were excluded included in our nearshore sampling analysis (i.e., blue crab and spot). If whole-subestuary mechanisms like TN or SAV are driving the results, offshore patterns for these species should mirror our nearshore results. Otter trawl samples were not included in the main local-or systemscale analyses because those analyses focused on nearshore areas.
Total nitrogen concentration (TN, μg L −1 ) was determined by analyzing water samples for nitrate by cadmium reduction and for ammonium plus organic nitrogen by Kjeldahl digestion; TN data were available from 27 subestuaries (collection dates 2010-2012) and concentrations ranged from 544 to 1050 μg L −1 . Within each subestuary, we averaged TN values from summer water samples collected at faunal sampling locations and from submerged aquatic vegetation sampling locations from a related study. Daytime dissolved oxygen was measured in subestuary channels offshore of a subset of sample sites and from sample sites with depth >2.5 using a YSI 600 QS (n = 101 sites and 12 subestuaries). Areal coverage of SAV, which typically occurs in waters <2 m deep, was estimated sensu Patrick et al. (2014) within each subestuary for each year from 2002 to 2012 and then averaged across years. This time period was chosen because it overlaps with the years in which fauna were sampled. Average SAV area was then converted to density by dividing average SAV area by subestuary area (i.e., km 2 of SAV per km 2 subestuary area).
We also conducted several a posteriori analyses to help address spatial confounding of land cover present in this study. In Chesapeake Bay, subestuaries with high amounts of hardened shoreline vs. high amounts of whole-catchment cropland tended to be on opposite shores (Fig. 1), and both stressors were positively related to distance from the ocean (R 2 = 0.59 for % watershed cropland and 0.63 for % hardened shoreline, second-order relationships). This spatial arrangement is not unusual because major cities and ports often develop at the upstream end of waters navigable by large shipping vessels (e.g., Baltimore, Washington D.C., London, Montreal), and agriculture grows near cities to satisfy the demand for food. Under such circumstances, identifying patterns between land cover and ocean-spawning species can be challenging because their abundances could decline with distance from the ocean. To address this issue, we included distance from the ocean as a predictor for all ocean-spawning species in the multivariable models used in the main analysis. We also compared the abundances of three ocean-spawning species (blue crab Callinectes sapidus, spot Leiostomus xanthurus, and Atlantic croaker) in 16 subestuary pairs that had comparable distances to the ocean (difference <20 km) but different amounts of % watershed cropland (difference > 10%). Pairwise data were analyzed using two-group Mann-Whitney U tests.

Local-Scale Effects of Shoreline Type
Abundances of 10 of 12 fish and crustacean species tested were strongly related to local shoreline type (beach, wetland [native vegetation], bulkhead, or riprap revetment; Fig. 3, Online Resources B), with specific effects dependent on functional species group (Table 1) and body size. Shoreline type explained 20.8 and 14.7% of the variation in faunal community structure within 0-3 and 0-16 m from shore, respectively (redundancy analysis, Online Resource C). Littoral-demersal species (small-bodied bottom-oriented species often found in shallow water) including mummichog (Fundulus heteroclitus), naked goby (Gobiosoma bosc), grass shrimp (Palaemonetes pugio), and striped killifish (Fundulus majalis) were more abundant adjacent to wetlands than hardened shoreline types at both distances from shore sampled (Fig. 3). Striped killifish and two planktivores (bay anchovy (Anchoa mitchilli) and Menidia spp.) were also more abundant at beaches within 0-16 m of shore than hardened shoreline types (only significant for striped killifish and Menidia spp.), while blue crab abundance and fish species richness were significantly lower at beaches than other shoreline types within 0-3 m from shore. In contrast, several benthivorepiscivores (larger-bodied bottom-oriented species) including spot, white perch (Morone Americana), and striped bass (Morone saxatilis), one open-water planktivore (bay anchovy), and Centrarchidae, which are either planktivores or benthivore-piscivores depending on life stage, were more abundant within 0-3 m of shore along hardened shorelines (especially bulkheads) than at wetlands or beaches (Fig. 3  upper panel). Within 0-16 m of shore, benthivore-piscivore species abundances were similar between hardened and wetland shorelines (Fig. 3 lower panel), suggesting some patterns within 0-3 m of shore reflect small shifts in habitat use. Notably, both the ANOVA analysis (here) and redundancy analysis (Online Resource C) showed similar species abundances on riprap and bulkhead habitat within 0-3 and 0-16 m from shore. See Online Resource B (Tables B1 and B2) for complete statistical output, including F and P values, of blocked ANOVAs examining local-scale effects of shoreline habitat on faunal abundances in waters within 0-3 and 0-16 m from shore. See Online Resource C (Fig. C1) for output of a redundancy analysis corroborating these results.
Local-scale water depth was significantly greater at bulkheads and ripraps compared to wetlands and beaches at 0 m and 3 m distance from shore (Online Resource B [ Table B3]; average of 0.51 and 0.26 m deeper, respectively). However, depth increased with distance from shore (Online Resource B, Fig. B1) and was similar across all four shoreline types at 16 m distance (shoreline type averages ranged from 0.99-1.07 m, p = 0.93). Two species associated with hardened shorelines in waters within 0-3 and 0-16 m from shore (spot and bay anchovy) and two larger-bodied species with comparable abundance at wetlands and hardened shoreline in waters within 0-3 and 0-16 m from shore (blue crab and Atlantic croaker) were actually most abundant along wetland shorelines in samples collected farther from shore (10 to 50 m; Fig. 4).

System-Scale Predictors of Fish and Crustacean Abundance
Watershed land cover, cumulative shoreline condition, and local shoreline type were prominent predictors of faunal abundance. Effects of these key predictors are summarized in Table 2; see Online Resource B (Tables B4 and B5) for complete output of system-scale models. The relative importance of each predictor depended on functional species group and dominant watershed land cover ( Table 2). Abundances of benthivore-piscivore and planktivore species were most influenced by system-scale predictors, especially % watershed cropland, the % of land within 100 m of a subestuary's shore comprised of wetlands (i.e., B% nearshore wetlands^), and the % of subestuarine shoreline comprised of hardened structures (i.e., B% hardened shoreline^) (Table 2; Fig. 5). These predictors occurred in the best models explaining species abundances after accounting for physical system characteristics  Local-scale associations of 12 fish and crustacean species with natural (wetland and beach) and hardened (bulkhead and riprap) shoreline types. Species are ordered to clump similar patterns. Abundances and species richness are reported in Z-score units after standardization with respect to gear type and year of capture; 0.0 on the y-axes are equal to average abundance for a species within each gear/year. Numbers above species codes on the x-axes correspond to statistically significant pairwise differences in habitat associations (See Online Resource B [Tables B1 and  B2] for complete statistical output). 1 = wetland > all others; 2 = wetland > bulkhead and riprap; 3 = wetland > beach; 4 = beach > all others; 5 = beach <all others, 6 = beach > bulkhead and riprap; 7 = bulkhead > all others; 8 = bulkhead > beach; 9 = bulkhead > wetland; 10 = riprap > beach; 11 = riprap > wetland (e.g., average subestuary water depth and surface area) and known environmental drivers (e.g., salinity and distance from the ocean, Online Resource B [Tables B4 and B5]). In contrast, abundances of littoral-demersal species were most influenced by local-scale shoreline type, possibly because of small bodies and relatively small home ranges (Lotrich 1975), and because of strong associations with shallow water that is less available near hardened shores.
The importance of watershed land use and shoreline hardening to faunal abundance differed between the Bcropland subset^of watersheds (excludes watersheds with >25% developed land cover) and the Bdeveloped land subset' of watersheds (excludes watersheds with >18% cropland; see Methods and Fig. 2). Relationships with % watershed cropland were significant for 6 of 15 species (4 negative, 2 positive) in the cropland subset (Table 2; Fig. 5), but predictors other than whole-watershed land use usually dominated species abundance models in the developed land subset. In particular, % hardened shoreline (a positive correlate of developed land use, Fig. 6) was significant for 10 of 15 species (9 negative, 1 positive) in the developed land subset (Table 2; Fig. 5). On average, subestuaries dominated by developed land (i.e., >25% developed land) had 45.2 ± 5.9% (±SE) hardened shoreline, compared to 16.7 ± 2.7% (±SE) for subestuaries dominated by cropland (>18% cropland) and 16.9 ± 2.2% (±SE) for 'control' subestuaries (i.e., subestuaries with ≤18% cropland and ≤25% developed land). Similarities in % hardened shoreline between high cropland and Bcontrolŝ ubestuaries further illustrate that % hardened shoreline in our study is largely a consequence of developed land use.
Faunal abundances were related to % nearshore wetlands in both land use subsets, with significant relationships for 13 of 15 species (10 positive, 3 negative; Table 2; Fig. 5). Although land cover was a direct predictor for fewer species, % nearshore wetlands was negatively correlated with developed land (R 2 = 0.17, p = 0.06) and cropland (negative threshold relationship). In addition, % nearshore wetlands and % hardened shoreline were inversely related (R 2 = 0.36, p < 0.0001) and thus are opposite measures of the same underlying variable (cumulative shoreline condition). Abundances of five species were affected by both cropland cover and cumulative shoreline condition. Three species (blue crab, spot, and Atlantic croaker) were negatively related to cropland and % hardened shoreline and positively related to % nearshore wetlands (Fig. 5). Similarly, mummichog was negatively related to cropland and % nearshore wetlands, but positively associated with wetlands at the local scale (see below for more on this apparent contradiction). By contrast, the species that supports   Table B1). Percent of shoreline occupied by emergent wetland vegetation was visually assessed at each site. Each data point represents the mean abundance from sites with a given percentage and error bars represent ±SE. All four relationships are significant at α = 0.05 the largest commercial fishery in Chesapeake Bay (Atlantic menhaden) was positively associated with cropland and % nearshore wetland. Multivariable models explained >50% of amongsubestuary variation for six of eight benthivore-piscivore species in both the cropland and developed land subsets (Table 2). Watershed land cover and cumulative shoreline condition accounted for much of this explanatory power, even when salinity or distance from the ocean were included as predictors.
The two main anthropogenic stressors identified by this study-% watershed cropland and % hardened shorelinemay have cumulative effects. To explore this possibility, we developed a composite stressor index for all 39 sampled subestuaries by adding watershed cropland and % hardened shoreline together using z-scored data, which have comparable units. We also developed an index of vegetated fish habitat-% nearshore wetlands + SAV-using the same method because these variables are negatively correlated with one another (Patrick et al. 2014). Subestuaries with aboveaverage stressor levels had below-average amounts of vegetated fish habitat and below-average abundances of three representative benthivore-piscivore species (blue crab, spot, and Atlantic croaker) (Fig. 7). Thus, watershed cropland and % hardened shoreline likely act simultaneously in negatively affecting bottom-oriented fishes and crustaceans, possibly because cropland and % hardened shoreline are both associated with loss of wetlands and SAV (Fig. 8). Relationships of three benthivore/piscivore species with the three critical system-scale predictors as identified by multi-variable mixed effects models. All three predictors were included in the best models or in alternate models (within 2.0 ΔAIC units) for these species (Table 2 and Online Resource B). Each point represents the average catch of the species within one subestuary, ±SE. Abundances are reported in Z-score units after standardization with respect to gear type and year of capture; 0.0 on the y-axes are equal to average abundance for a species within each gear/year. Blue crab (Callinectes sapidus), spot (Leiostomus xanthurus), and Micropogonias undulates (Atlantic croaker, abbreviated to BCroaker^) were chosen to illustrate these patterns because of their ecological and economic importance in Atlantic coastal waters and because of their functional similarity as benthivore/piscivores. Patterns with % watershed cropland included subestuaries from the cropland subset, patterns with % hardened shoreline included subestuaries from the developed land subset, and % nearshore wetlands patterns included all subestuaries
Several possible mechanisms for patterns with watershed cropland were evaluated a posteriori, including total nitrogen, dissolved oxygen, SAV and % nearshore wetlands. The relative importance of each mechanism varied by species (Table 3). Positive relationships between cropland and two planktivores (Atlantic menhaden and Centrarchidae [although Centrarchids may also be benthivore-piscivores]; Table 2) appeared to be driven by TN, which explained more amongsubestuary variation than cropland, SAV, or % nearshore wetlands (Table 3). By contrast, negative relationships between cropland and three benthivore-piscivores (Table 2) appeared to be driven primarily by SAV density, which explained more among-subestuary variation than watershed cropland and TN for Atlantic croaker and spot (Table 3). However, models containing SAV, TN, and their interaction in place for cropland explained the most among-subestuary variation for Atlantic croaker and spot, while SAV and TN explained comparable amounts of among-subestuary variation in blue crab abundance (Table 3). Taken together, these analyses suggest both SAV and TN may be contributing to patterns between benthivore-piscivore abundance and watershed cropland, with SAV playing a stronger role. % nearshore wetlands was also important for all three benthivore-piscivores, especially Atlantic croaker and also explained the most amongsubestuary variation in the abundance of the littoraldemersal species (mummichog).
Distance from the ocean was spatially confounded with land cover (i.e., positive relationships between distance and both % cropland and % developed land), but several a posteriori analyses indicated that % watershed cropland and cumulative shoreline condition had real effects on oceanspawning species that went beyond relationships with distance from the ocean. For example, three ocean-spawning species (blue crab, spot, and Atlantic croaker) had significantly greater abundance in subestuaries with lower % watershed cropland when comparing pairs of subestuaries that were comparable distances to the ocean (Fig. 9). In addition, the oceanspawning planktivorous Atlantic menhaden was positively related to watershed cropland (Table 2) even though % watershed cropland increased with distance from the ocean. Finally, our multivariable approach identified % watershed cropland and % hardened shoreline as important model components in addition to distance from the ocean for several ocean- . Cropland data are plotted on the cropland subset and developed land data on the developed land subset. TN data are averaged from water samples collected during summer during faunal and SAV sampling, and error bars are ±SE. DO data are daytime measurements from subestuary channels offshore of samples sites and measurements from sample sites with depth >2.5 m. Percent hardened shoreline data are derived from the Comprehensive Coastal Inventory Program and do not have variance (hence, no error bars). p = 0.002, 0.005, and 0.23 for the relationships between cropland and TN, DO, and % hardened shoreline, respectively. p = 0.40, 0.18, and <0.0001 for the relationships between developed land and TN, DO, and % hardened shoreline, respectively spawning species (Table 2) suggesting effects beyond those expected with distance. Although distance from the ocean certainly has an effect on these species, the weight of the evidence suggests that % watershed cropland and cumulative shoreline condition both have real effects on mobile fish and crustaceans.

Discussion
The results of this study provide the most comprehensive picture to date of how shoreline hardening and land use affect mobile fish and crustacean abundances. We were able to identify potentially elusive relationships because we studied a large number of sites and systems in Chesapeake Bay and the Delaware Coastal Bays and included both juvenile and adult life stages. Our results indicate both shoreline hardening and high levels of agricultural land use are likely to affect the abundance of a number of common fish and mobile crustacean species, including several that support commercial and recreational fisheries, in nearshore estuarine habitat. We add to a growing body of evidence that shoreline hardening strongly affects marine faunal communities at local scales. Moreover, combined with results of recent, large-scale research in Puget Sound (Dethier et al. 2016), which differs dramatically from Chesapeake Bay in temperature, species composition, and its extensive rocky shorelines, it is clear shoreline hardening can have cumulative impacts on physical, geomorphic, and biological elements of nearshore ecosystems.

Local-Scale Effects of Shoreline Type
Hardened shorelines generally had lower abundances of small-bodied fishes and greater abundances of larger-bodied fishes at the local scale, especially within 3 m from shore. Faunal abundances were also similar between the two hardened shoreline types (i.e., bulkhead and riprap revetment) with only one significant difference (striped bass within 3 m of shore were significantly higher at bulkhead). This is notable because riprap revetment is generally considered to have lower impact than bulkhead, a notion that is reflected in the environmental policy of some states (e.g., North Carolina statute 15A NCAC 7H .0208 [b][7] (2012), which states Bwhere possible, sloping rip-rap, gabions, or vegetation shall be used rather than bulkheads^). Although this may be true for some species or processes, our study suggests riprap revetment is R 2 = 0.42 R 2 = 0.33 R 2 = 0.31 R 2 = 0.47 7 Fig. 7 Relationships of vegetated fish habitat and the abundances of three benthivore/piscivore species with a composite index of cropland plus % hardened shoreline. Each point represents the average catch of the species within one subestuary, ±SE (n = 39 subestuaries). Blue crab (Callinectes sapidus), spot (Leiostomus xanthurus), and Micropogonias undulates (Atlantic croaker) were chosen to illustrate these patterns because of their ecological and economic importance in Atlantic coastal waters and because of their functional similarity as benthivore/piscivores.
Abundances are reported in Z-score units after standardization with respect to gear type and year of capture; average values (i.e., z-score = 0) are noted with dotted lines. Both the composite stressor index (cropland + % hardened shoreline) and vegetated fish habitat (SAV + % nearshore wetland) were developed by adding together z-scored values, which have comparable units. Model fit (linear vs. exponential decay) was determined by maximizing R 2 ; non-linear models were only used if they improved R 2 by at least 0.05 not better than bulkhead for mobile fish and crustaceans. This is a key finding relevant to shoreline policies.
Results of our comprehensive study substantiate previous smaller-scale research programs that found reduced densities of some small-bodied fishes and crustaceans along hardened shorelines (Peterson et al. 2000;Rozas et al. 2007;Peterson & Lowe 2009;Balouskus & Targett 2016). Small-bodied fauna may favor natural shorelines at local scales because they have more resources than hardened shorelines. Compared to other shoreline types, wetlands provide larger amounts of allochthonous carbon to adjacent nearshore waters (Weinstein et al. 2005), support greater densities of infaunal and epifaunal invertebrates (King et al. 2005;Seitz et al. 2006), and are associated with more subtidal benthic structure (e.g., woody debris, shellfish beds) that may serve as habitat (Bilkovic & Roggero 2008). Similarly, natural beaches have more organic material and greater density and diversity of benthic invertebrates than armored beaches (Sobocinski et al. 2010), and beach alteration can negatively impact nursery areas (Peterson & Bishop 2005). Munsch et al. (2015b) also found reduced epibenthic prey in the water column at hardened shorelines when compared to beaches, which in turn affected the diets of some fish. Greater wave reflection at hardened shorelines may also decrease habitat quality for small-bodied fishes by increasing turbulence, scouring/damaging the benthos, and reducing water clarity through fine sediment resuspension (Pope 1997;Miles et al. 2001).
The local-scale shoreline type associations observed in this study may also reflect added effects of hardening on predatorprey dynamics. Hardened shorelines foster greater local-scale water depth in near shore areas (Online Resource B, Fig. B1), reducing shallow water refuge for small-bodied species and R 2 = 0.27 R 2 = 0.37 Fig. 8 Negative relationships between vegetated fish habitat (a composite index of SAV and nearshore tidal wetlands) and both % watershed cropland cover (top, p = 0.006) and % hardened shoreline (bottom, p = 0.004). Cropland data are plotted on the cropland subset, and % hardened shoreline data on the developed land subset. SAV and nearshore wetlands were presented as a composite index, calculated by adding z-scored values of SAV and % nearshore wetland together, because both forms of vegetation provide fishes with critical refuge habitat and because SAV and % nearshore wetland are negatively correlated with one another (Patrick et al. 2014) affording large-bodied species greater access to nearshore areas (Hines & Ruiz 1995). Shoreline type-specific differences in water depth at 0 and 3 m from shore, compared with similar water depths at 16 m from shore, may explain the observed differences in faunal responses at different distances from shore. The positive effect of hardened shoreline on largebodied or open-water species abundances is more likely a preference for deeper water than for hardened shoreline structure, supporting other studies (Toft et al. 2007). Small-bodied species likely face increased predation risk along hardened shorelines at the local scale, given known positive relationships between nearshore predation and water depth ). Our empirical evidence (Figs. 3 and B1) supports the hypothesis that human alteration of nearshore bathymetry via shoreline hardening negatively affects localized abundances of species dependent on shallow water as refuge habitat (Hines & Ruiz 1995).

Cumulative, System-Scale Effects of Shoreline Condition
Our study offers compelling evidence that cumulative shoreline condition is a key predictor of faunal abundance at the system scale. When effects were detected (10 of 15 and 13 of 15 species in the cropland and developed land use subsets, respectively), most relationships were negative regarding % hardened shoreline and positive regarding % nearshore wetlands. Seven of 15 species were positively related to % nearshore wetland in the cropland subset, and 9 of 15 (and fish species richness) were positively related to % nearshore wetland in the developed land subset. Conversely, 5 of 15 and 9 of 15 species were negatively related to % hardened shoreline in the cropland and developed land subsets, respectively. In total, the abundance of 13 of 15 species, as well as fish species richness, were related to cumulative shoreline condition (either % hardened shoreline or % nearshore wetland) in at least one of the two land use subsets (Table 2). This evidence strongly suggests the loss of natural shorelines and the construction of hardened shorelines affects mobile fish and crustaceans at the system scale and is consistent with recent findings that shoreline hardening has cumulative effects on geomorphic processes, such as beach grain size and slope, that likely affect habitat quality for fish and invertebrates (Dethier et al. 2016). Moreover, our findings provide empirical evidence that substantiates concerns about the negative cumulative impacts of shoreline hardening (Peterson & Lowe 2009) and the extensive global destruction of natural shorelines like tidal wetlands (Lotze et al. 2006;Barbier et al. 2011) on fauna.
The effect of cumulative shoreline condition was relatively consistent across functional groups. System-scale abundances of 5 of 8 benthivore-piscivores were negatively related to % hardened shoreline in at least one of the two landuse subsets vs. 2 of 4 planktivores and 2 of 4 littoral demersals. Only Table 3 Testing species abundance relationships to potential mechanistic predictors for species in which % watershed cropland was a significant predictor by (A) evaluating the % of among-subestuary variation explained by substituting mechanistic predictors for % cropland in multivariable models and (B) evaluation the % of among-subestuary variation explained by models including only one predictor We explicitly test two mechanistic predictors associated with % watershed cropland (submerged aquatic vegetation, hereafter BSAV^, and total nitrogen, hereafter BTN^) when substituted for % watershed cropland in species abundance models. In addition, we report the % of among-subestuary variation explained by models that contained only the random effect of subestuary of capture and one of five fixed-effect predictors: % watershed cropland, SAV, TN, % nearshore wetland, and distance from the ocean. Subestuaries with >25% of the watershed comprised of developed land were excluded. Average total nitrogen data were available for 18 such subestuaries, and we constrained our analysis to this subset even though SAV data were available from more subestuaries so that each predictor was assessed on the same dataset. Model predictors were those from the best model in Online Resource B (Table B4) except as noted: 1 One of the alternate models for Atlantic croaker abundance identified in Online Resource B (Table B4) had the lowest AIC score on this subset of data and, thus, was used for means of comparison (included % watershed cropland and % nearshore wetlands as predictors). 2 One of the alternate models for spot abundance identified in Online Resource B (Table B4) had the lowest AIC score on this subset of data and, thus, was used for means of comparison (included % watershed cropland, % nearshore forest, % hardened shoreline, and distance from the ocean as predictors). 3 % cropland was included in alternate model of Atlantic menhaden abundance (the alternate model included only % watershed as a predictor) and, thus, this species is tested here. *a single high-influence subestuary (in this case, a subestuary with extremely high menhaden abundance compared to all other subestuaries) was removed from % variation calculations Centrarchids, which are more freshwater oriented throughout their lifecycle than any other species in our study, increased in abundance with increasing shoreline hardening in both land use subsets; this pattern may be an artifact of developed land use usually occurring near the head of each subestuary in the Chesapeake Bay, where salinities are typically low (Peterson & Meador 1994). Hogchoker, the only other species that responded positively to increasing shoreline hardening, did so only in the cropland subset, where it also responded positively to increasing wetland shoreline; hogchoker responded negatively to hardened shoreline in the developed land subset.
Similarly, system-scale abundances of 6 of 8 benthivorepiscivores were positively related to % nearshore wetlands in at least one of the two land use subsets, vs. 3 of 4 planktivores and 1 of 4 littoral-demersals. Replacing wetlands with hardened shorelines likely has cumulative system-scale effects on benthivore-piscivores by reducing the abundance of benthic invertebrate prey (Seitz et al.2006) and by simplifying subtidal benthic structure (Bilkovic & Roggero 2008). Wetland habitats also benefit some planktivorous species, including those positively associated with % nearshore wetlands in this study, by providing spawning and nursery areas (Love & May 2007;Balouskus & Targett 2012) and by providing allochthonous nutrients and organic matter to support secondary production (Weinstein et al. 2005). Wetland habitats are also undoubtedly important for littoral-demersal species, but this appeared to primarily be a local-scale effect: all four littoral-demersals in our study were positively associated with wetland habitat at the local scale (Fig. 3) despite limited positive relationships with % nearshore wetland at the system scale.
The system-scale effects of cumulative shoreline condition were generally consistent with effects of local-scale shoreline type. For example, striped killifish abundance was positively associated with beach shoreline at the local scale (Fig. 3) and the % of shoreline comprised of beach at the system scale ( Table 2). In addition, several benthivore-piscivores (blue crab, spot, and Atlantic croaker) and one planktivore (bay anchovy) were positively associated with wetland shoreline at the local scale in samples collected 10 to 50 m offshore (Fig. 4) and at the sytem scale (Fig. 5). However, some littoral-demersal species such as mummichog and striped killifish had greater local-scale abundance at wetland shorelines (Fig. 3), but were negatively correlated with % nearshore wetlands at the system scale (Table 2). This contradiction could reflect greater dispersal of littoral-demersal species within more extensive wetlands; perhaps larger wetlands allowed such species to remain on wetland surfaces or in tidal creeks instead of occupying the nearshore zone.
Watershed land cover was less universally important than shoreline condition in predicting system-scale faunal abundances and, in this study, was only a direct predictor of faunal abundance in watersheds in the cropland subset. This difference may relate to the mechanisms by which shoreline condition and land cover affect fishes and crustaceans (Fig. 10).

Blue Crab
Spot Atlantic Croaker Fig. 9 Pairwise comparisons of blue crab (top), spot (middle), and Atlantic croaker (bottom) abundance in subestuaries with comparable distances to the ocean (difference <20 km) but different amounts of % watershed cropland (difference ≥10%). Based on relationships between these species and % watershed cropland (Fig. 5), Blow^cropland subestuaries had ≤20% watershed cropland while Bhigh^cropland subestuaries had >20%. There is a visible decline in the abundance of all three species as distance from the ocean increases (left to right). Abundances in the subestuary with less cropland cover were ≥ abundances in the subestuary with more cropland cover for 13 of 16 pairs for blue crab (p = 0.03), 12 of 16 pairs for spot (p = 0.05), and 14 of 16 pairs for Atlantic croaker (p = 0.007). Pairwise data analyzed using twogroup Mann-Whitney U Tests. Please see Online Resource A (Table A1) for details on subestuaries used in this comparison Wetlands provide habitat for many nearshore species (Beck et al. 2001) and are biological filters of nutrients, contaminants, and sediments (Vernberg 1993) independent of land cover; thus, it is no surprise that shoreline condition was important in both the cropland and developed land subsets. By contrast, the effects of cropland and developed land on adjacent aquatic habitat varied in the subestuaries we analyzed. Cropland cover was associated with reduced water quality as well as SAV and nearshore wetlands, while developed land was more strongly correlated with shoreline hardening (Fig. 6). Shoreline hardening may be a symptom of developed land use (Gittman et al. 2015) but may also mask other developed land effects. For example, total nitrogen and dissolved oxygen have been previously linked with impervious surfaces and developed land use by other studies (Foley et al. 2005;Rabalais et al. 2009;Uphoff et al., 2011), but not in the sites analyzed in the present study. Notably, shoreline hardening was negatively associated with SAV and wetlands and may cause SAV loss by increasing wave reflection, which can directly affect SAV through scour, and indirectly by reducing water clarity through sediment suspension (Pope 1997;Miles et al. 2001).

System-Scale Effects and Potential Mechanisms of % Watershed Cropland
Our cross-system analysis indicates that the abundances of several key species are related to % watershed cropland. Bottom-oriented benthivore-piscivores such as blue crab and Atlantic croaker had lower abundances in systems with high cropland cover, while planktivorous species such as Atlantic menhaden had higher abundances in cropland dominated systems. These relationships may be driven by several potential mechanisms (Fig. 10). Cropland is strongly tied to non-point nutrient pollution (i.e., eutrophication) in coastal systems (Jordan et al. 1997), which can increase plankton production and lead to hypoxia (Kemp et al. 2005;Diaz & Rosenberg 2008) and can reduce coverage and density of SAV by decreasing water clarity and promoting growth of epiphytic algae (Fisher et al. 2006;Patrick et al. 2014). We found evidence for this in our study, where increasing % watershed cropland was linked with nutrient enrichment and reduced dissolved oxygen as well as scarcity of SAV and nearshore wetlands.
Our results provide strong empirical evidence for the assertion that watersheds dominated by croplands affect mobile fish and crustaceans at the system-scale through multiple mechanisms related to eutrophication. A posteriori analyses provided evidence that SAV density and % nearshore wetlands were likely the strongest mechanisms behind negative relationships between % watershed cropland and benthivorepiscivore abundance. These findings are consistent with earlier studies demonstrating that SAV benefits benthivorepiscivore species including blue crab, Atlantic croaker, and spot by providing food resources and structural habitat/ predation refuge (Laney 1997;Heck et al. 2003). Total nitrogen also appeared to contribute to negative benthivore-piscivore/cropland patterns to a lesser extent, but appeared to be the primary driver of positive associations between % watershed cropland cover and the abundances of two planktivore species. Nitrogen concentrations per se would not directly affect fish, crabs or grass shrimp. Rather, positive Fig. 10 Conceptual diagram of mechanistic relationships between land cover, shoreline hardening, and the abundance of estuarine fish and crustacean from three generalized functional groups. Positive and negative signs correspond with the direction of each relationship. Large, blocked arrows represent main effects found by this study.
Narrow, dashed arrows represent mechanistic relationships sometimes shown to be important by other studies, but not directly observed by this study. Other than these two distinctions, arrow size is not representative of effect size. The three generalized functional species groups are circled relationships between planktivore abundance, cropland cover, and TN are consistent with studies showing that planktivores respond positively to high nutrient loads because of greater amounts of plankton (Caddy 1993;de Leiva Moreno et al., 2000). Elevated TN levels also likely contribute to negative relationships between benthivore-piscivores and cropland by stimulating planktonic and epiphytic algal blooms (Fisher et al. 2006) and thereby reducing water clarity, which in turn reduces SAV density and impedes foraging efficiency for visual predators (Henley et al. 2000). This argument is substantiated by negative interactions between SAV and TN observed in models for all species whose abundance was related to cropland (significant for spot and Atlantic croaker). Blue crab, spot, and Atlantic croaker (as well as many other species) annually migrate into and out of the studied subestuaries and, thus, abundance patterns may be related to avoidance of subestuaries with poor habitat and water quality (Breitburg et al. 2009). High TN loads might also negatively affect benthivore-piscivores by exacerbating deepwater (Kemp et al. 2005;Breitburg et al. 2009) or shallow diel-cycling (Tyler et al. 2009) hypoxia, which can reduce benthic food supply via losses of sessile prey (Diaz & Rosenberg 2008); alter fish behavior (Brady & Targett 2013); or impair the reproductive capacity of estuarine fishes (Cheek et al. 2009). This hypothesis is supported by lower daytime bottom dissolved oxygen (DO) in high-cropland subestuaries (Fig. 6, middle panel). However, both cropland and developed land had similar relationships with daytime bottom DO (Fig. 6, middle panel). This suggests that rather than hypoxia, another aspect of elevated TN (such as algal blooms and reduced SAV) drove relationships between faunal abundances and cropland.
Patterns with land use may be more likely to emerge at the system (subestuary) or watershed scale than the local scale. A recent Chesapeake Bay study noted no effect of the percent of developed and agricultural land use within 500 m of sample sites on abundances of several fish and crab species in seagrass beds (Blake et al. 2014). By contrast, our study identified several significant patterns with cropland use at the watershed scale. Otter trawl samples of mid-channel bottom habitat from a subset of subestuaries demonstrated that blue crab and spot from offshore areas also showed negative relationships to % watershed cropland cover (R 2 = 0.30 for blue crab and 0.15 for spot), supporting the hypothesis that observed cropland effects are related to whole-subestuary mechanisms like elevated TN or reduced SAV rather than localized effects in nearshore areas.

Study Limitations
Several limitations must be acknowledged when interpreting this study. First, although we argued that the patterns in species-specific abundances shown by our study are related to beneficial or detrimental qualities of surrounding habitat, abundance is not necessarily an indicator of habitat value. For example, fish abundances may also be related to transient occupation of certain areas due to life histories and/or spatial arrangement of habitat types. Additional research into why species select for or survive better in certain habitats would be beneficial. Second, our study did not separately examine adult and juvenile life stages, and assumed that catch vulnerability of adults and juveniles was equal. Although this assumption was probably true for many species in this study, different life stages of large-bodied animals (e.g., striped bass, white perch, Atlantic croaker) may have had different gear vulnerabilities or ecological responses. Size-class specific utilization of shoreline type will be explored in a future study. Third, our study only covered the range of conditions that were available and could be safely sampled by the methods used in this study. Although the subestuaries included in our study are representative of land use and shoreline type gradients prevalent in Chesapeake Bay and elsewhere, there are no truly pristine subestuaries in Chesapeake Bay because of tidal exchange with the mainstem. We also avoided sampling some of the most polluted subestuaries due to health concerns related to discharge of untreated sewage. Fourth, there were insufficient data to thoroughly explore possible roles of sediment or contaminant loads, which are also associated with cropland (Foley et al. 2005). However, sediment or contaminant loads would be unlikely to positively affect abundance of planktivorous Atlantic menhaden and juvenile Centrarchids as observed in our study. Fifth, by limiting our study to only the most commonly occurring species, we were able to evaluate patterns across broad gradients but may have missed details relevant to assemblage composition as well as functional, taxonomic, and beta diversity; these topics are the subject of concurrent studies that include both common and uncommon species. Sixth, our study compiled data from two dominant time periods (2002-2003 and 2009-2012). In order to focus on relationships between species abundances, shoreline hardening, and land use, we deliberately excluded interannual effects on species abundances through our z-scoring standardization procedure, and there are likely interannual effects present in the unstandardized data that may be informative for studies investigating temporal patterns. Finally, although our use of z-scored abundances follows well-known methods for meta-analyses combining data from different sampling protocols (e.g., Felker-Quinn et al. 2013;Shade et al. 2013), their use makes abundance values less tangible. However, differences in z-score abundance represent meaningful differences in the sampled abundance of each species in our analysis, which was limited to the most common species and included over 600,000 individuals. For example, blue crab density ranged from 0 to 0.7 m −2 , and spot from 0 to 5.4 m −3 , in one of datasets included in the analysis.

Conclusions
Until the present study, solid linkages between fish and crustacean abundance, watershed land cover, and cumulative shoreline condition have been elusive, albeit long suspected (e.g., Kemp et al. 2005;Toft et al. 2007;Peterson & Lowe 2009;Dethier et al. 2016). Our study provides clear evidence that shoreline hardening affects local fish and crustacean assemblages, and that cumulative shoreline hardening negatively affects fauna at the system scale. We expect that effects on fish and crustaceans will worsen as shoreline alteration intensifies in coastal areas due to burgeoning human populations (Crossett et al. 2004;Brody et al. 2008) and sea level rise (Rahmstorf 2007;Peterson & Lowe 2009). Despite concern regarding shoreline habitat alteration, public perception has been slow to change, and management recommendations and regulations have been inconsistent (Peterson & Lowe 2009). Controlling shoreline hardening will likely benefit ecologically and economically important species, but will be challenging as sea level rise threatens housing and economically important land, requiring action from resource managers and policymakers that can navigate difficult decisions at local and regional scales. For example, property owners often decide to armor shorelines to prevent damage from armoring on neighboring property (Scyphers et al. 2015); this would suggest that coordinated local-scale efforts would be needed to prevent or reverse shoreline hardening at large spatial scales.
Our empirical cross-system study also helps clarify the ongoing debate regarding the effects of eutrophication on populations of mobile fish and crustaceans by demonstrating different effects on planktivores vs. benthivorepiscivores, as proposed by others (Caddy 1993;de Leiva Moreno et al., 2000). We show that croplands negatively affect the abundance and/or distribution of several bottom-oriented marine species in shallow nearshore waters. Although only 4 of 15 species were negatively related to % watershed cropland, these species support highly valuable commercial and recreational fisheries (blue crab, spot, and Atlantic croaker) or are a key prey item in nearshore foodwebs (mummichog). A posteriori analysis demonstrating that SAV and/or TN explained more variation than did cropland cover in the abundance of these species suggests that nutrient loading was probably responsible for the observed patterns. This would also suggest that measures to reduce nutrient loads may benefit bottomoriented fishes and crabs with high ecological/economic value and may also have negative effects on fisheries of planktivorous species.
Finally, we show that the nearshore abundances of most species are strongly related to factors other than land cover. Measures to reduce nutrient loadings would have many benefits (Rabalais et al. 2009), but should not be viewed as a panacea for protecting fauna. For example, species whose abundances were negatively related to cropland cover in our study were also strongly related to cumulative shoreline condition (positive relationships with % nearshore wetland and negative relationships with % hardened shoreline). Moreover, relationships with land cover were found only with cropland, while cumulative shoreline condition was an important predictor of the abundance of most species across all subestuaries. Shoreline condition and land cover likely have synergistic effects given that natural shorelines buffer against runoff of nutrients and other pollutants (Vernberg 1993). However, land cover and eutrophication have been targeted by recovery and conservation efforts (USEPA 2010), while shoreline hardening has received less attention (Peterson & Lowe 2009). Nutrient load reduction initiatives and policies should continue, but fauna in estuaries and other coastal systems would also benefit from increased efforts to limit construction of hardened structure, pursue less harmful coastal defense strategies such as living shorelines (e.g., Gittman et al. 2016), and conserve natural shorelines.