Human Influence at the Coast: Upland and Shoreline Stressors Affect Coastal Macrofauna and Are Mediated by Salinity

Anthropogenic stressors can affect subtidal communities within the land-water interface. Increasing anthropogenic activities, including upland and shoreline development, threaten ecologically important species in these habitats. In this study, we examined the consequences of anthropogenic stressors on benthic macrofaunal communities in 14 subestuaries of Chesapeake Bay. We investigated how subestuary upland use (forested, agricultural, developed land) and shoreline development (riprap and bulkhead compared to marsh and beach) affected density, biomass, and diversity of benthic infauna. Upland and shoreline development were parameters included in the most plausible models among a candidate set compared using corrected Akaike’s Information Criterion. For benthic macrofauna, density tended to be lower in subestuaries with developed or mixed compared to forested or agricultural upland use. Benthic biomass was significantly lower in subestuaries with developed compared to forested upland use, and biomass declined exponentially with proportion of near-shore developed land. Benthic density did not differ significantly among natural marsh, beach, and riprap habitats, but tended to be lower adjacent to bulkhead shorelines. Including all subestuaries, there were no differences in diversity by shoreline type. In low salinities, benthic Shannon (H′) diversity tended to be higher adjacent to natural marshes compared to the other habitats, and lower adjacent to bulkheads, but the pattern was reversed in high salinities. Sediment characteristics varied by shoreline type and contributed to differences in benthic community structure. Given the changes in the infaunal community with anthropogenic stressors, subestuary upland and shoreline development should be minimized to increase benthic production and subsequent trophic transfer within the food web.


Introduction
Coastal ecosystems are threatened by anthropogenic stressors as human populations flock to coastal areas in record numbers (Halpern et al. 2007;Airoldi and Beck 2007). These coastal habitats are highly productive and serve important roles as feeding grounds, nursery areas, spawning areas, and corridors for migration of ecologically and commercially important marine species (Beck et al. 2001;Seitz et al. 2014); therefore, the loss or degradation of these habitats can have dramatic effects on ecosystem productivity. Among coastal areas, estuarine habitats and wetlands are especially productive, serving multiple ecosystem functions, with their value expected to increase in the future (Costanza et al. 1997(Costanza et al. , 2014. In response to sea-level rise and coastal erosion, estuarine shorelines are being hardened at alarming rates, with losses to wetlands and coastal habitats (Gittman et al. 2015). Multiple stressors, including the combination of upland and shoreline development, could have synergistic effects on estuarine fauna (Crain et al. 2008). The combination of these factors on estuarine fauna has rarely been examined (King et al. 2005;Li et al. 2007;Patrick et al. 2014), and never, to our knowledge, for the infaunal, deep-dwelling benthic community.

Upland Use
Construction of human infrastructure has resulted in increased watershed development and runoff (Jantz et al. 2005), which

Communicated by Carles Ibanez Marti
Electronic supplementary material The online version of this article (https://doi.org/10.1007/s12237-017-0347-6) contains supplementary material, which is available to authorized users.
can add sediments, nutrients, pesticides, and contaminants to water (Jordan et al. 1997;Paul and Meyer 2001;Gregg et al. 2015). This can negatively impact benthic invertebrates (Hale et al. 2004;King et al. 2005), fish (Sanger et al. 2004), and waterbird communities (DeLuca et al. 2008;Prosser et al. 2017). Pollutants can reduce abundance, diversity, and trophic complexity of benthic species in favor of small, short-lived, opportunistic species, such as deposit-feeding polychaetes (Pearson and Rosenberg 1978;Warwick and Clarke 1993;Inglis and Kross 2000). Examining the combined effects of shoreline and upland development will help managers incorporate trade-offs between anthropogenic use and impacts on estuarine fauna into management and conservation decisions.

Shoreline Alteration
The extensive use of bulkheads (vertical seawalls) and riprap (rocky revetments) to armor shorelines and the construction of jetties, docks, piers, and marinas all result in the replacement of natural habitat with artificial structures, with modifications to the dynamics of these critical systems (Gittman et al. 2016b). Shoreline armoring has increased in developed countries including the Netherlands, Japan, and the USA; for example 14% of the US coastline has been armored, with some subestuaries having over 50% of shorelines developed (Gittman et al. 2015). In other areas (e.g., France, Spain, Italy), 45% of the coastal zone has been developed (Bulleri and Chapman 2010), and erosion protection through shoreline development is a problem (Jiménez et al. 2016;Santana-Cordero et al. 2016;Harik et al. 2017).
Shoreline hardening is of particular interest in Chesapeake Bay. The dearth of information on the ecological effects of shoreline structures (Weinstein and Kreeger 2000) can limit managers' understanding of habitat degradation, reducing their ability to make environmentally sound decisions. In Chesapeake Bay, it remains unknown whether there are thresholds of shoreline and upland development that, if exceeded, lead to loss of ecosystem services, but thresholds of development have been identified in other systems (Dethier et al. 2016). With shoreline development, there may be some species that are Bwinners^(e.g., opportunistic species) and some that are Blosers^(e.g., sensitive species) (Weisberg et al. 1997); thus, overall diversity may not necessarily change with shoreline development. Therefore, our study aims to quantitatively estimate the effects of upland and shoreline development on the shallow, subtidal benthos of Chesapeake Bay.

Salinity
Estuarine salinity gradients structure flora and fauna, particularly in benthic communities, and they can mediate effects of anthropogenic stressors. For example, differences in salinity result in differential responses of seagrass to shoreline development in Chesapeake Bay (Patrick et al. 2016). Moreover, effects of physical or chemical stressors, such as hypoxia, have differential effects on benthos depending on salinity regime (King et al. 2005;Seitz et al. 2009). Predators may be more abundant in high-salinity areas, and coupling of the benthos with adjacent habitats may be greater in low-salinity areas, leading to differential effects on benthos in different salinity regimes. For example, blue crabs and fish are more abundant in high-salinity subestuaries (King et al. 2005), particularly near marsh habitats (Kornis et al. 2017), and their feeding may reduce benthic biomass (Menge and Sutherland 1987). Low-salinity (≤ 15 psu), upper-Bay subestuaries are typically smaller, shallower water bodies than high-salinity subestuaries, leading to closer coupling of the benthos with adjacent habitats that would be linearly related to benthic exchange with the water column (Gerritsen et al. 1994). Therefore, responses of benthic communities to upland and shoreline development may be expected to differ by salinity regime.

Faunal Responses
Shoreline development and upland use affect benthic communities and predators. Benthic communities are effective indicators of ecological condition and harbingers of ecological stress (Widdicombe and Spicer 2008). Abundance and diversity of subtidal, infaunal benthic invertebrates are higher adjacent to natural marsh than bulkhead shorelines, intermediate at riprap shorelines, and predator density and diversity tends to be higher adjacent to natural marsh shorelines . Negative effects of hardened shorelines can be compounded when there is extensive development in the surrounding landscape (Seitz and Lawless 2008).
Further, there can be wider, ecosystem-level impacts of development. A Chesapeake Bay-wide trawl survey identified benthic invertebrates as the predominant source of carbon in the diets of fishes (Buchheister and Latour 2015). Given the key role of benthic invertebrates in the food web, ecologically and economically important species of invertebrates and finfish may be stressed by both habitat loss and prey reduction at the land-water interface. Several studies report negative effects of altered shorelines on predators in adjacent waters Hendon et al. 2000;Carroll 2003;Peterson and Lowe 2009;Kornis et al. 2017), and one review summarizes the negative effects on estuarine fish (Munsch et al. 2017), but these studies did not concurrently examine infaunal benthos. biomass, and diversity and thereby exacerbate the effects of shoreline hardening, and (2) shoreline development reduces local density, biomass, and diversity of benthic infauna. We compared density, diversity, and biomass of benthic macrofauna adjacent to four shoreline types across a range of salinities and land uses in replicate subestuaries of Chesapeake Bay. We also tested the generality of previous findings on the effects of shoreline alterations through a large-scale, multisubestuary empirical study. The experimental design included (i) two salinity regions-high salinity (generally polyhaline to meso-polyhaline), or > 15 psu, and low salinity (generally low-mesohaline), or ≤ 15 psu-to capture a range of conditions and account for salinity-driven differences among regions, (ii) subestuaries of three differing predominant land usage patterns (forested, agricultural, or developed), and (iii) four shoreline types per subestuary (natural Spartina marsh, sandy beach, riprap, and bulkhead [= seawall]) and thereby assess the impact of multiple stressors on benthic communities (Table 1).

Study Locations
The Chesapeake Bay is 300-km long, and its shoreline consists of over 100 subestuaries (Li et al. 2007;Weller and Baker 2014). Each of these subestuaries is unique, and their watersheds have differing proportions of forested, agricultural, and developed upland use (Li et al. 2007;Patrick et al. 2014). Urban development pressure tends to be highest in the upper Bay (King et al. 2005) due to the large cities in that area. The eastern shore of the Bay is heavily developed with agriculture, including crops and chicken farms, leading to runoff of sediments, nitrogen, and phosphorous (Jordan et al. 1997;Prasad et al. 2014).
We investigated 14 subestuaries throughout the Chesapeake Bay from 2010 to 2013 (Fig. 1). The subestuaries were selected based on their primary upland use, forested, agricultural, developed or mixed, as characterized by Li et al. (2007) and Patrick et al. (2014), and on the amount of hardened shoreline and shoreline condition (VIMS-CCRM: http://ccrm.vims.edu/gisdatabases.html). Dominant watershed land cover was obtained from the 2006 National Land Cover Dataset, which was derived from Landsat 7 satellite remote sensing imagery with 30-m resolution (Fry et al. 2011). The watersheds were each classified (by Li et al. 2007 andPatrick et al. 2014) into a category of dominant land cover based on the following: B(1) forested (≥ 60% forest and forested wetland), (2) developed (≥ 50% developed land), (3) agricultural (≥ 40% cropland), (4) mixed-developed (15-50% developed land), (5) mixed-agricultural (20-40% cropland), and (6) mixed-undisturbed (watersheds which did not fit into categories 1-5)^ (Patrick et al. 2014). We used the term Bmixed^for category 6. The subestuaries included developed subestuaries (Stony Creek, Magothy River, Mill Creek, and Poquoson), agricultural subestuaries (Miles River, Harris Creek, Onancock Creek, Occohannock Creek), forested subestuaries (Monroe Bay, Corrotoman River, East River, Severn River, and Catlett Islands), and a mixed subestuary (Yeocomico Creek) (Fig. 1). Two to seven different subestuaries were sampled annually between late June and early August in 2010 to 2013 ( Table 2). The percentage of watershed use was not available for Catlett Islands, so it was interpolated using the mean calculated from the upland-use values from each of two nearby creeks, Poropotank Bay and Queen's Creek. We chose subestuaries within the mesohaline to polyhaline regimes (5-30 psu) so that the benthic communities were not composed of oligohaline species.

Survey Methods
In each subestuary, four shoreline types were sampled: marsh, beach, riprap revetment, and bulkhead. Four to six replicates were sampled at each shoreline in each subestuary. As in other surveys of this type, there may have been some spatial confounding because shoreline types were clustered within subestuaries. In 2010, four randomly selected replicates were collected at each shoreline type resulting in 16 samples in both East River and Occohannock Creek. From 2011 to 2013, six replicate samples were collected at each of the four shoreline types resulting in 24 samples in each subestuary. However, in three subestuaries (Harris Creek, Mill Creek, and Catlett Islands), samples were collected at only 3, 1, and 0 beaches, respectively, due to the limited number of beach shorelines. Where available, we used areas that had > 30 m of the k, the number of parameters in each model, which includes 1 for variance. If a β is located in a column, then that variable was included in the model. Excluding the Yeocomico, models were run for all subestuaries, lowsalinity subestuaries (Sal ≤ 15 psu), and high-salinity subestuaries (Sal > 15 psu) particular shoreline type, and we randomly selected from among multiple sites of a given shoreline type. Five to six replicates per treatment are sufficient to detect differences in densities of infauna among shoreline types . Infaunal organisms were collected five meters from shore using a benthic suction sampler designed to capture deepdwelling macrofauna (Eggleston et al. 1992), and this distance from shore has been sufficient to show significant effects of shoreline type on benthic organisms in past research . We made sure to take 2-4 samples of each shoreline type each on rising and falling tides (mean tidal range of 1 m). We used a 0.11-m 2 PVC cylinder inserted to a depth of 30 cm in the sediment, evacuated the contents of the cylinder, and sieved it through a 3-mm mesh bag. This sampling targets large macrofauna, such as bivalves, that are deep dwelling, sparsely distributed, and are also important prey items for blue crabs and epibenthic fish (Seitz et al. 2003). The contents of the mesh bag were frozen until sorting.
In the laboratory, each sample was sorted thoroughly and double-checked with a second sorting. Infauna were identified to the lowest taxon possible (usually species except some polychaetes, e.g., Capitellidae and Spionidae), enumerated, and stored in 70% ethanol. Organisms were dried for 48 h at 70°C and then combusted in a muffle furnace for 4 h at 550°C to obtain ash-free dry weights (AFDW). For each sample, bulk weights were obtained for most taxa (e.g., polychaetes, crustaceans, anemones), while bivalve species were weighed separately.
Our research examined direct and indirect effects of multiple stressors and focused on benthic communities as well as ecologically important individual benthic species. Direct effects of shoreline development and upland use include  (circle, agricultural; square, developed; triangle, forested; octagon, mixed) changes in macrofaunal species composition. Indirect effects include changes in sediment grain size and composition, lowered dissolved oxygen, reduced submerged aquatic vegetation (SAV) abundance, and the spread of invasive saltmarsh plants in ways that secondarily influence macrofauna. The ecologically important species in our study included some long-lived, pollution-sensitive (= sensitive) taxa (sensu the Chesapeake Bay Benthic Monitoring Program's categorization: Weisberg et al. 1997), such as the clam Limecola balthica. This is an important sessile prey species that cannot alter its distributions in response to stressors, unlike fish and mobile invertebrates. Sensitive taxa included the polychaetes, Clymenlla torquata and Glycera americana, and the clams, Limecola balthica, Rangia cuneata, and Tagelus plebeius, which are indicative of a mature community (Weisberg et al. 1997;Llansó et al. 2002). We also examined responses of pollution-indicative (= tolerant) taxa, including the small polychaetes Leitoscoloplos spp. and Eteone heteropoda (Dauer 1993), and the clam Mulinia lateralis (Weisberg et al. 1997;Llansó et al. 2002), which are opportunistic, weedy species that are most common in degraded habitats (Pearson and Rosenberg 1978;Dimitriou et al. 2015). We compared tolerant taxa and sensitive taxa among shoreline types.
Water quality was measured with a calibrated YSI Pro-Plus Multiparameter at each site for temperature (°C), salinity, and dissolved oxygen (mg 1 −1 ). Additionally, two sediment samples were collected with a 4-cm 2 syringe to a depth of 5 cm. One core was collected for total organic carbon and nitrogen (TOC/TN) content, and an Exeter CE440 elemental analyzer was used to quantify carbon, hydrogen, and nitrogen (CHN) content. The second core was used for grain-size analysis using a standard wet sieving and pipetting technique (Plumb 1981).

Statistical Analyses
The effects of environmental variables (shoreline type, upland use, salinity, and sediment type) on infaunal density, biomass, and Shannon diversity were examined in an information theory framework using general linear models. We used simple linear regression to examine relationships between some variables to establish independence. TOC/TN were positively related to each other and with sediment type (see Lawless and Seitz 2014) and therefore could not be included as independent variables in the analyses. Density estimates were obtained by standardizing the raw infaunal abundance per replicate by the surface area of the PVC cylinder, and we calculated a mean and standard deviation from the six replicates to obtain density, biomass, or diversity per shoreline type. Shannon (H′) diversity and richness (S) were calculated using PRIMER v 7.
Explanatory mathematical general linear models were created based on multiple working hypotheses regarding influential variables. Akaike's Information Criterion (AIC) was used to evaluate the log likelihood of each explanatory model within a candidate set of models while accounting for the number of parameters of each model (Anderson 2008). Likelihood was estimated from general linear models and raw data were transformed if they deviated from a normal distribution. To account for the potential bias associated with a small sample size, the corrected AIC (AICc) was used; the correction factor approaches zero as the sample size increases (Anderson 2008). Parameter estimates and standard errors were Table 2 Year sampled and mean values (± SE) of water quality and sediment of each subestuary (with the upland use type indicated in parentheses)

Subestuary
Year generated via multimodel inference across all models containing a particular parameter (Anderson 2008). We proposed a total of eight candidate models (Table 1) for each of the three response variables. In addition to using AICc to explore the full dataset, we also examined model selection for the data split into low-(salinity ≤ 15 psu) and high-(salinity > 15 psu) salinity subestuaries. The high-salinity subestuaries included the East, Occohannock, Poquoson, Severn, Catlett, and Onancock rivers, and low-salinity subestuaries included Corrotoman River, Harris Creek, Magothy River, Miles Creek, Monroe Bay, Stony Creek, and Mill Creek. Our division agreed with classification of these subestuaries as high salinity (polyhaline) or low salinity from long-term Chesapeake Bay monitoring (Li et al. 2007). The Yeocomico subestuary was excluded from AIC analyses, as it was the only subestuary with mixed upland use. While diversity and biomass followed a normal distribution, density exhibited a right-skewed distribution and was log-transformed. For AIC analyses, samples containing zero individuals (3 instances out of 284 samples) were removed. Analyses were run in the open-source statistical package R (R Development Core Team 2016).
We also used general linear models in some cases to more closely examine differences among levels of categorical upland or shoreline factors separately for three different response variables (benthic density, biomass, diversity), with α ≤ 0.05 as the cutoff for statistical significance. Fisher post-hoc multiple comparison tests were used to determine differences between levels of each factor for significant general linear models. We used non-linear least squares regression to examine mean community biomass versus proportion of land developed within 250 m of the shoreline across all subestuaries (Fig. 2), since non-linear relationships were hypothesized based on previous work on macrofaunal responses to developed land use (Bilkovic et al. 2006).
Distance-based permutational multivariate analysis of variance (PERMANOVA; McArdle and Anderson 2001;Anderson 2001) was used to test for differences in benthic communities by shoreline type and upland use. Bray-Curtis similarities were calculated on square-root transformed species abundance and biomass matrices (to normalize the data) and permutated 9999 times. The design for the analysis consisted of three factors: shoreline (Sh), four levels, fixed; upland use (Up), three levels, fixed, and river (Ri), nested in upland (Up) 13 levels, fixed. We used the Type III sum of squares within PERMANOVA to determine significance ) because we had an unbalanced sampling design with different numbers of subestuaries in each upland type.
Non-metric multidimensional scaling (nMDS) ordination plots of centroids were used to summarize patterns of benthic assemblages by the factors upland use (forested, agricultural, and developed) and shoreline (marsh, beach, riprap, and bulkhead). The centroid was determined as the center point of all samples for a certain shoreline within an upland type in multidimensional space. Matrices of the distances among centroids were derived from Bray-Curtis similarity matrices of square-root transformed infaunal community data. The centroid ordination plots were derived from the distances among centroid matrices generated from Bray-Curtis matrices of square-root transformed abundance and biomass data. All PERMANOVA, nMDS, and distance between centroid species contributions were calculated in the PRIMER v 7 PERMANOVA+ add-on package Clarke et al. 2014).

Physical Characteristics
Temperature and dissolved oxygen varied little among subestuaries, but salinity varied substantially among subestuaries (Table 2). Sediment % sand + gravel was generally high, as 13 of the 14 subestuaries had sediment > 80.0% sand + gravel. Overall by shoreline type, beaches had the highest percentage of sand + gravel (94.60%), followed by Fig. 2 Upland-use percentage 250 m from shore for experimental subestuaries. The Bagricultural^category includes cultivated crops, such as corn and soy, and includes grasses and hay riprap (88.46%), bulkhead (86.94%), and lastly marsh (80.37%). Many forested and agricultural subestuaries (e.g., Miles River, Monroe Bay, Severn River, Corrotoman River) tended to have high TOC (Table 2), whereas developed or mixed subestuaries (e.g., Magothy River, Yeocomico River) had some of the lowest TOC. Both total organic carbon (TOC) and total nitrogen (TN) were two times higher at marsh than at bulkhead and riprap and four times higher than at beaches. TOC and TN were significantly and tightly related (R 2 = 0.94), and both were inversely related to % sand + gravel (R 2 = 0.46 and 0.30, respectively).

Upland Use
We collected 37 benthic taxa in the 3-mm suction samples throughout the Bay. Polychaetes were the dominant taxon by richness (21 taxa) followed by bivalves (11 species). Thirteen taxa contributed to 90% of the abundance. The nereid polychaete Alitta succinea contributed most to abundance (28.92%), followed by the clam Limecola balthica (13.71%), spionid polychaetes (10.33%), and the clam Ameritella mitchelli (9.60%). By upland use, similar taxa were dominant across the subestuary upland-use types, but the percent contributions of each taxon to the community varied, with 4-8 taxa encompassing over 90% of total abundance (Table 3 ). Agricultural subestuaries were dominated > 70% by four taxa, A. mitchelli, Alitta succinea, L. balthica, and capitellid polychaetes, and these four cumulatively contributed to the same percentage of abundance as did the polychaete A. succinea alone in forested subestuaries. Developed and mixed-developed subestuaries were also dominated by A. succinea, and secondarily by Ameritella mitchelli and Rangia cuneata in developed subestuaries, and spionids in the mixed-developed subestuary.
Comparing density among differing upland-use types across all subestuaries, there were no significant differences but some notable trends. Density tended to be higher in forested and agricultural subestuaries than in developed or mixed subestuaries (Fig. 3a). In low salinity, agricultural subestuaries tended to have the highest densities, forested and developed subestuaries were intermediate, and the mixed subestuary tended to be lower (Fig. 3b). In high salinity, densities tended to be higher in forested, intermediate in agricultural, and lowest in developed subestuaries (Fig. 3c).
There were significant differences in biomass among upland-use types. Biomass was significantly higher in forested and agricultural subestuaries than in developed or mixed subestuaries (Fig. 3a, d)   high salinity, biomass was not statistically different across upland-use types, with a tendency for highest biomass in forested subestuaries (Fig. 3f).
Overall, richness and diversity were relatively low, with highest mean values of 6 and 1.4, respectively. There were some significant differences in richness and diversity among subestuary upland-use types. Richness among all subestuaries was highest in agricultural compared to forested, developed, and mixed subestuaries (Fig. 4a) (general linear model and Fisher test p < 0.0001). In low-salinity regions, richness was highest in agricultural and developed subestuaries (Fig. 4b) (general linear model and Fisher test p < 0.0001) while in high-salinity regions, richness was not statistically different across upland-use types (Fig. 4c). In all subestuaries, lowand high-salinity regions, Shannon (H′) diversity was higher in agricultural and developed subestuaries than in forested subestuaries ( Fig. 4d-f) (general linear model and Fisher test p < 0.0001, p < 0.0001, and p = 0.003, respectively). Across subestuaries, there was an exponential decline in infaunal biomass with % development in the zone 250 m from the shore (Fig. 5). In the subestuaries with < 20% upland development, mean biomass ranged from a high of 10.5 g AFDW m −2 to a low of 1.5 g AFDW m −2 . In subestuaries with ≥ 20% development, biomass remained low and never reached higher than 3.5 g AFDW m −2 .
Community assemblages differed by upland use. Centroid nMDS ordination plots of abundance assemblages by shoreline type in each upland type clustered tightly by upland use and did not overlap other upland-use groups ( Supplementary  Fig. A1a). Centroid nMDS ordination plots of biomass assemblages clustered by upland use, and shoreline types for forested and developed uplands did not overlap other upland-use Fig. 3 Upland-use effects on mean (± SE) 3-mm infaunal a density of all subestuaries, b density of subestuaries ≤ 15 psu, c density of subestuaries > 15 psu, d biomass of all subestuaries, e biomass of subestuaries ≤ 15 psu, and f biomass of subestuaries > 15 psu. Upland type: For, forested; Ag, agricultural; Dev, developed; and Mix, mixed. Means are from replicates among all shoreline types within a subestuary. Small letters above bars denote significant differences determined by a Fisher post-hoc multiple comparison test at α = 0.05. NA indicates data not available groups. Centroids of shoreline types within forested upland use clustered the tightest while centroids of shoreline types within agricultural upland use had greater dispersion. The centroids of agricultural marsh and beach biomass were closer to the forested shoreline types than to agricultural beach and riprap shoreline types. (Supplementary Fig. A1b). PERMANOVA revealed significant effects for BUpland Use,^BShoreline,^and BRiver^(nested within BUpland Use^) for abundance assemblages (Table A1). A significant interaction between BShoreline^and BRiver^(nested within BUpland Use^) indicated that shoreline effects varied between rivers and by upland use. For biomass assemblages, PERMANOVA revealed significant effects for BUpland Use,^BShoreline,^and BRiver^(nested within BUpland Use^) (Table A1).

Shoreline Development
Across all subestuaries, there were no significant differences in density by shoreline type, though some trends were evident (Fig. 6a-c). Among all subestuaries, density tended to be lowest at bulkhead shorelines (Fig. 6a). In low-salinity subestuaries, density tended to be higher at beaches than all other shoreline types and lowest at bulkheads (Fig. 6b). In high-salinity subestuaries, density did not differ among shoreline types and no consistent patterns were evident ( Fig. 6c; general linear model: p = 0.378).
Biomass did not differ significantly by shoreline type, and trends were mixed (Fig. 6). Among all subestuaries, biomass was equivalent among bulkhead, riprap, and beach habitats (Fig. 6d). In the low-salinity subestuaries, biomass tended to Fig. 4 Upland-use effects on mean (± SE) 3-mm infaunal a richness of all subestuaries, b richness of subestuaries ≤ 15 psu, c richness of subestuaries > 15 psu, d diversity of all subestuaries, e diversity of subestuaries ≤ 15 psu, and f diversity of subestuaries > 15 psu. Upland type: For, forested; Ag, agricultural; Dev, developed; and Mix, mixed. Small letters above bars denote significant differences determined by a Fisher post-hoc multiple comparison test at α = 0.05. NA indicates data not available be higher at beaches than all other shoreline types (Fig. 6e). In the high-salinity subestuaries, biomass did not differ among shoreline types but trends were opposite of low-salinity subestuaries ( Fig. 6f; general linear model: p = 0.184).
Across all subestuaries, richness and Shannon diversity did not differ; they tended to have opposite patterns by shoreline type depending on whether they were high-salinity or lowsalinity subestuaries (Fig. 7). In the low-salinity subestuaries (salinity ≤ 15 psu), richness and Shannon diversity tended to be higher at marshes and beaches than developed shoreline types (Fig. 7b, e). In the high-salinity subestuaries (salinity > 15 psu), richness and diversity were lowest at marshes (general linear model p = 0.018 and 0.005, respectively).

Taxon-Specific Responses to Shoreline Development
Individual taxon patterns by shoreline type (Table 4) were different depending on whether the taxon was long-lived and sensitive or short-lived and opportunistic (tolerant) (pollution sensitivity defined by Weisberg et al. 1997). Limecola balthica biomass was highest at marshes and beaches; biomass at marshes was nearly two times higher than at riprap and three times higher than at bulkheads (Fig. 8a) (general linear model: p = 0.041). Ameritella mitchelli (a sensitive clam) biomass tended to be higher at natural than hardened shorelines (Fig. 8b). Mulinia lateralis (a tolerant clam) biomass tended to be lowest at marshes while highest at beaches followed by riprap and bulkhead (Fig. 8c). Rangia cuneata (a tolerant clam) biomass tended to be highest at hardened shorelines (riprap and bulkhead) and lowest at marshes; bulkheads had five times higher biomass than marshes (Fig. 8d). Further, some taxa were only present at natural shorelines (e.g., Glycera americana) or were absent at bulkheads (e.g., Cirratulidae and Edotea triloba).
The percentage of tolerant taxa tended to be highest at bulkheads while the percentage of sensitive taxa tended to be lowest at bulkheads (Fig. 9a, b). Two tolerant polycheate taxa, E. heteropoda and Leitoscoloplos spp., tended to have highest densities at bulkheads, followed by riprap, beaches, and marshes (Fig. 9c, d).

Combined Development Influences
Based on the AIC model weights, the model that only contained upland use was typically the top or one of the top models for all response variables (Table A2). For density, top models contained upland use, shoreline, and/or sediment type, and the global model was supported in the low-salinity grouping. Multimodel inference indicated that macrofaunal density and biomass varied among subestuaries with different upland use, with developed subestuaries having significantly lower density and biomass than agricultural or forested subestuaries (Table A3). For biomass, the model with upland use was consistently supported in all three subestuary groupings, and it was the best model for total and low-salinity subestuaries. For low-and high-salinity subestuaries, supported models also included shoreline, and for total subestuaries, the null model was also supported. For diversity, the global model was the best model for total subestuaries; however, for low salinities, the model with upland use was the best model, and for high salinities, the shoreline plus upland use plus sediment model was supported. Multimodel inference indicated that macrofaunal diversity varied among subestuaries with different upland use, with agricultural and developed subestuaries having significantly higher diversity than forested subestuaries (Table A3), though the number of taxa differed by only a few (Fig. 4a).

Multiple Anthropogenic Stressors
The multiple stressors of upland and shoreline development influenced benthic infaunal communities throughout Chesapeake Bay (Tables A1 & A2). Anthropogenic impacts can act in concert to negatively impact benthic communities. Our results agree with previous studies in the lower Chesapeake Bay where the effects of shoreline development on macrofauna depended on co-occurring stressors, such as shoreline hardening and upland use (Seitz and Lawless 2008;Davis et al. 2008;Peterson and Lowe 2009).

Upland Use
Across all subestuaries examined, benthic communities adjacent to upland watersheds that were forested or agricultural tended to have higher density and had significantly higher biomass than those adjacent to developed subestuaries or the Fig. 5 Percentage of development within 250 m of shoreline versus mean infaunal community biomass per subestuary (g ash-free dry weight; AFDW) with exponential decline curve estuary with mixed development (Figs. 3d and 5). For density, this pattern was driven by significant differences in highsalinity subestuaries; however, for biomass, both high-and low-salinity subestuaries showed similar patterns. Reduced biomass may be an indirect effect of the increased inflow of nutrients, toxicants, and sediments due to increased runoff over hardened surfaces (Jordan et al. 1997;Gregg et al. 2015), which can negatively affect some benthic organisms, especially species with low tolerance to stressful conditions. Our work concurs with previous studies in Chesapeake Bay that demonstrated the importance of both upland use and salinity for a few key species, namely Callinectes sapidus, Limecola (formerly Macoma) balthica, and Ameritella (formerly Macoma) mitchelli (King et al. 2005), but we expand this concept to the importance for the large, infaunal benthic community. In our study, infaunal biomass declined exponentially with upland development within 250 m of the shoreline, declining dramatically through 20% development, suggesting food availability decreased for higher trophic levels as upland development increased. Moreover, in other studies, landscape-level effects masked shoreline effects (Seitz and Lawless 2008;Lawless and Seitz 2014); therefore, upland effects are paramount. Densities of tolerant taxa increased with near-shore upland development, as the abundances of small, opportunistic taxa increased. The detrimental effects of pollutants associated with urban uplands likely contributed to increased variability with increased upland development, and this increased variability has been shown previously in stressed benthic communities (Warwick and Clarke 1993 Resource control of benthic communities may be affecting the overall trends of higher benthic density and biomass with reduced upland development. Benthic food availability can affect distributions of benthic species (Diaz and Schaffner 1990), and food for deposit-feeders may be of higher quality adjacent to both forested and agricultural watersheds where natural, high-carbon, allochthonous material runs off, as compared to developed watersheds with less allochthonous carbon input (Dauer et al. 1992;Rodil et al. 2008). Sedimentary food availability (organic carbon and nitrogen) for benthic organisms, which increases with lower % sand, was higher in many of the rivers with forested and agricultural upland use that we sampled (e.g., Miles River, Monroe Bay, Severn River). This potentially contributed to higher benthic densities of deposit-feeding species (e.g., A. mitchelli, L. balthica, and Alitta succinea; Lovall et al. 2016) in those locations.

Shoreline Development
In terms of shoreline development, bulkhead habitats tended to have reduced benthic density compared to natural habitats across all subestuaries. This is likely because marshes act as efficient filters of runoff from the upland (Howes et al. 1996;Roman et al. 2000), but bulkheads sever that land-water buffer, allowing excess toxicants to enter the water and deter benthic organisms. Several studies suggest that shoreline development may strongly affect macrofauna in nearby subtidal shallow zones (Tourtellotte and Dauer 1983;Weis et Davis et al. 2008;Peterson and Lowe 2009). This investigation encompassing the entire Chesapeake Bay agrees with earlier work showing detrimental effects of hardened shorelines within individual subestuaries Bilkovic and Roggero 2008;Gittman et al. 2016a;Lovall et al. 2016), but our large-scale study across a diverse set of subestuaries advances our knowledge about where, within a large estuary, effects are most prominent, and which species are most severely affected. Specifically, effects of shoreline development on density and richness were most prominent in low-salinity subestuaries. The pattern of increased density with natural shorelines was driven by large, long-lived sensitive taxa, such as the clams L. balthica and A. mitchelli (Weisberg et al. 1997;Long et al. 2014), which can deposit feed and take advantage of detrital material delivered from marshes (Kamermans 1994). The response of organisms to riprap structures was somewhat intermediate between bulkhead and natural habitats, as the unconsolidated rock structures possibly provided some improved habitat that was not provided by bulkhead habitats. Some taxa were less sensitive to pollutants that easily runoff from impervious surfaces of developed shorelines (Jordan et al. 1997;Paul and Meyer 2001;Gregg et al. 2015); tolerant taxa (Weisberg et al. 1997), such as Eteone heteropoda and Leitoscoloplos spp., had highest densities adjacent to bulkhead and riprap shorelines, which lack the filtering capabilities of wetland buffers (Howes et al. 1996;Roman et al. 2000). These taxa are typically considered Bweedy^and indicative of deteriorated habitat conditions (Weisberg et al. 1997), suggesting that bulkhead and riprap shorelines reduce functionality of benthic habitats.
The trend of biomass by shoreline was somewhat counter to our hypothesis, possibly because of sediment effects on some abundant high-biomass, suspension-feeding bivalves (e.g., Tagelus plebeius, Rangia cuneata, and Mya arenaria; Table 4). These species were deterred in muddier sediments of the natural marsh habitats where a bivalve's feeding apparatus can get clogged with fine sediment (Steele-Petrović 1975), an indirect effect of shoreline development. There may have also been some effect of increased predator abundance at natural marsh habitats Bilkovic and Roggero 2008;Kornis et al. 2017), which could have resulted in lower benthic biomass there during our mid-summer sampling when predation peaks (Moody 2001). Upland, shoreline, and salinity were important in predicting benthic biomass, possibly driven by large bivalves (e.g., the introduced species Rangia cuneata) present particularly in developed, low-salinity subestuaries (Table 3) and Tagelus plebeius and Mya arenaria in high-salinity developed subestuaries and at bulkhead shorelines (Figs. 6f and 8d). Upland effects were paramount, as benthic community structure differed with upland type (Fig.  A1a), and accounting for upland differences, community structure differed by shoreline within rivers. Biomass of indicators of high habitat quality, e.g., Limecola balthica (Pearson and Rosenberg 1978;Long et al. 2014), was greatest adjacent to natural marsh habitats, but these taxa were only common in forested and agricultural subestuaries. Larger, long-lived indicator species are indicative of good ecosystem functioning, as chronic disturbance can reduce their abundance and decrease ecosystem functioning (Pearson and Rosenberg 1978). The tendency for natural marsh habitats to have higher benthic richness and Shannon (H′) diversity than bulkhead and riprap habitats in low-salinity subestuaries, and the reverse trend in high-salinity subestuaries, may have been driven by differences in predation and coupling with uplands. Predators such as blue crabs and fish are more abundant in high-salinity subestuaries (King et al. 2005), particularly near marsh habitats, as determined by a companion study of fish abundance in relation to shoreline development in Chesapeake Bay (Kornis et al. 2017). Excessive feeding by these predators may contribute to the loss of benthic biomass (Menge and Sutherland 1987). Prominent effects of shoreline development in lower-salinity (≤ 15 psu), upper-Bay subestuaries are possibly because they are more closely coupled with influences from adjacent habitats, as water depths are shallower with greater exchange of water with the benthos (Gerritsen et al. 1994). However, only a small number of taxa changed overall, and the species-specific changes were more noteworthy. This is the first study to explicitly examine both upland and shoreline development effects on benthic communities in relation to salinity, as the effects of multiple stressors on estuarine fauna have rarely been examined (King et al. 2005;Li et al. 2007;Patrick et al. 2014), and such studies have not examined the entire infaunal, deep-dwelling benthic community. Moreover, the increase in upland and shoreline development throughout the world (Gittman et al. 2015) has increased the demand for studies examining the interacting effects of these anthropogenic stressors. Ecological stress in benthic communities was detected in subestuaries with high percentages of urban development and in response to shoreline development, particularly for sensitive species (e.g., Limecola balthica, Fig. 8a). Our extensive spatial coverage throughout Chesapeake Bay and examination of multiple stressors was instrumental in demonstrating where development had the greatest effects (in low salinity, upper Chesapeake Bay), and that upland use was most influential on benthic communities, though shoreline development was also important. Further management efforts need to be aimed at reducing both upland and shoreline development to maintain the production of benthic communities that are important in the estuarine food web.