Individual, population, and ecosystem effects of hypoxia on a dominant benthic bivalve in Chesapeake Bay

Hypoxia is an environmental stressor that affects abundance, biomass, diversity, and ecosystem function of benthic assemblages worldwide, yet its collective impact at individual, population, and ecosystem levels has rarely been investigated. We examined the effects of hypoxia on the biomass-dominant clam, Macoma balthica, in the York and Rappahannock Rivers (Chesapeake Bay, USA). We (1) surveyed the M. balthica populations in both rivers in 2003 and 2004, (2) determined the effects of low dissolved oxygen (DO) onM. balthica fecundity in a laboratory experiment, and (3) employed a predator-exclusion field experiment to establish the effects of hypoxia and prey density on predation upon M. balthica. The resultant data were used to parameterize a matrix model, which was analyzed to define potential effects of hypoxia at the population level. In both rivers, hypoxia decreased individual clam growth and caused local extinction of populations. Hypoxia reduced egg production of M. balthica by 40% and increased protein investment per egg. In the predatorexclusion field experiment, hypoxia magnified predation rates threefold and altered the functional response of predators to M. balthica from a stabilizing type III functional response to a destabilizing type II functional response. In a density-independent matrix model, hypoxia resulted in coupled source–sink metapopulation dynamics, with hypoxic areas acting as blackhole sinks. Increases in the spatial and temporal extent of hypoxia caused the populations to decline toward extinction. In a second model that incorporated density dependence, under mild hypoxic conditions trophic transfer from M. balthica to predators increased, but at increased spatial or temporal extent of hypoxia trophic transfer decreased. The major decline in trophic transfer to predators under severe hypoxia resulted from diversion of M. balthica biomass into the microbial loop. Our model predicted that there are multiple stable states for M. balthica populations (high and very low densities), such that the saddle point (threshold at which the population switches from one state to the other) increased and resilience decreased with the spatial extent of hypoxia. We underscore how effects of a stressor at the individual level can combine to have substantial population and ecosystem-level effects.


Hypoxia and its ecological consequences
Anthropogenic eutrophication of estuarine systems has substantial negative effects on coastal ecosystems (Kemp et al. 2005). One important consequence of eutrophication is the development of hypoxia, or low dissolved oxygen (DO), in bottom waters Rosenberg 1995, 2008). Hypoxia usually develops when excessive algal growth, fueled by nutrient inputs, results in the flux of large amounts of labile organic material to bottom waters. If respiration depletes the oxygen faster than it can be replenished, this leads to the development of hypoxia, (DO , 2.0 mg O 2 /L; hereafter mg/L) or anoxia (DO ¼ 0.0 mg/L; Officer et al. 1984, Seliger et al. 1984, Rabalais et al. 2010.
Hypoxia can be a severe stressor in aquatic systems (Diaz and Rosenberg 2008, Breitburg et al. 2009, Levin et al. 2009), particularly for infaunal benthic species that are unable to migrate from the affected area. Severe hypoxia can kill benthic species and diminish abundance, biomass, diversity, and ecosystem function (e.g., Dauer et al. 1992, Powers et al. 2005, Montagna and Froeschke 2009. Also of interest are its sublethal effects at the individual level. Hypoxia can reduce individual feeding rates (Sagasti et al. 2001), growth (e.g., Nilsson 1999, Wiklund andSundelin 2001), reproduction (e.g., Condon et al. 2001, Wu andOr 2005), and antipredator responses (Wang et al. 2010). In addition, behavioral responses, such as reduced burial depth and extension of siphons or palps, are common in benthic infauna and can enhance susceptibility to predation (Taylor and Eggleston 2000, De Goeij et al. 2001, Seitz et al. 2003b, Montagna and Ritter 2006, Long et al. 2008b, Lee et al. 2012. Although much work has been conducted on the effects of hypoxia at the individual level, few studies explicitly consider how these effects interact over multiple scales to affect the population or ecosystem, though this is an important step in understanding largescale patterns (Diaz and Rosenberg 1995, 2008, Altieri and Witman 2006, Breitburg et al. 2009). In this study, we investigated the effects of hypoxia on the growth, reproduction, and survival of the Baltic macoma, Macoma balthica, in the York and Rappahannock Rivers, Chesapeake Bay, USA. We used these results to explore population-level effects by constructing a matrix model that examines how hypoxia-induced source-sink population dynamics affect population stability. We considered the effects that this will have on higher and lower trophic levels by performing a caging experiment to determine effects of hypoxia and prey density on the rate of predation. Finally, we incorporated the findings into a densitydependent model and explored the net results on ecosystem function and services.

Study area and species
Our studies were performed in the mesohaline region of the York and Rappahannock Rivers (Chesapeake Bay, Virginia, USA) in unstructured soft-bottom habitats (Fig. 1). Both rivers experience episodic hypoxia during the summer Rosenberg 1995, 2008). Typically, hypoxia in the Rappahannock is more severe than in the York River, with lower DO levels, a greater areal extent, and a longer duration (Kuo and Neilson 1987).
M. balthica, a small (shell length , 40 mm) thinshelled clam, is the biomass-dominant species in unvegetated soft-bottom habitats of Chesapeake Bay, where it can comprise over 85% of the biomass (Holland et al. 1977). With an LT 50 (lethal time for 50% of the population) of 8.3-18.0 days under nearanoxic conditions in the laboratory, it is highly tolerant of hypoxia (Henriksson 1969, de Zwaan et al. 2001. Under hypoxic conditions, this clam migrates to a shallower burial depth in the sediment and extends its siphons (Brafield 1963, Seitz et al. 2003b, Long et al. 2008b, which makes it more vulnerable to predation (De Goeij et al. 2001) because its primary defense against predation is deep burial Comtois 1985, Seitz et al. 2001). M. balthica is a key link in the food web by contributing to the flow of biomass from primary producers to upper trophic levels, such as the commercially important blue crab Callinectes sapidus (Baird and Ulanowicz 1989). In Chesapeake Bay, M. balthica comprises ;50% of the adult C. sapidus diet (Hines et al. 1990, Lipcius et al. 2007. Piscine predators, such as the Atlantic croaker (Micropogonias undulatus), spot (Leiostomus xanthurus), and hog-choker (Trinectes maculatus), are generally unable to reach whole clams due to deep clam burial, but will feed extensively on M. balthica siphons (Pihl et al. 1992, Derrick and Kennedy 1997, Powers et al. 2005, which can induce M. balthica to migrate vertically to a shallower burial depth (De Goeij et al. 2001).

Source-sink population dynamics
Source-sink theory recognizes that species occupy interconnected habitats differing in quality and carrying capacity within a heterogeneous landscape. These habitats can be separated broadly into sources and sinks. In sources, births exceed deaths and emigration exceeds immigration; population growth is positive and excess individuals emigrate. In sink habitats, the death rate exceeds the birth rate, and the populations are only sustained by immigration from source habitats (Pulliam 1988). In this way, a species can extend its range into areas where, without immigration, the population would suffer local extinction.
In a recent review of the evidence for source-sink dynamics in marine and estuarine species, there were only a few species with strong evidence for source-sink dynamics, but there was extensive circumstantial evidence for its widespread existence and a consequent need for further examination of the phenomenon (Lipcius and Ralph 2011). The most prevalent mechanisms underlying source-sink dynamics included variation in habitat quality (whether due to natural or anthropogenic causes), dispersal, predation, and fishery exploitation. It was further concluded that investigation of source-sink dynamics is critically needed to promote effective conservation and restoration of marine and estuarine species (Lipcius and Ralph 2011), and thus the emphasis upon source-sink dynamics in our study.
It may seem counterintuitive that an organism would live in a suboptimal (i.e., sink) habitat, but this is likely to be common in terrestrial (Pulliam and Danielson 1991) and marine (Lipcius and Ralph 2011) ecosystems. High predator abundance in source habitats can make it worthwhile for prey to move temporarily into sink habitats (Schreiber et al. 2006). Sink habitats can be populated due to competition forcing individuals out of the source habitat, or through passive dispersal from the source into the sink as might occur for benthic species with pelagic larval stages (Dias 1996). When the latter occurs, sink habitats can reduce the equilibrium population density of the source populations (Pulliam and Danielson 1991).
Many sedentary marine invertebrates, such as M. balthica, are susceptible to source-sink dynamics due to dispersal processes that transport larvae to unfavorable habitats (Lipcius and Ralph 2011); substrate types, salinity variations, current structure, and food availability are examples of factors that could cause an area to be a sink. Even though M. balthica larvae are capable of discerning the quality of habitats to some degree (Beukema and de Vlas 1989, Armonies 1992, Seitz et al. 2001, settlers can only detect current conditions and not future adverse conditions. As M. balthica recruits in the fall and spring, and hypoxia occurs in the summer, settlers are unlikely to distinguish future hypoxic areas.
Thus, we hypothesized that hypoxic areas of Chesapeake Bay tributaries act as population sinks that may negatively influence populations in normoxic source habitats.

Predation and density dependence
The functional response (FR) of a predator relates the rate of predation to the density of prey. Three general types of functional response are: type I, or linear and density independent; type II, or hyperbolic, which is inversely density dependent; and type III, or sigmoid, which is density dependent (Holling 1959, Hassell et al. 1977. Type II FRs are common where the prey species uses armor to increase handling time and reduce the rate of predation (Seitz et al. 2001). As the density of prey increases, the proportional mortality decreases, resulting in a high-density refuge from predation. When prey avoid predators by hiding or using camouflage, thus increasing the searching time and decreasing the predation rate, type III FRs are common (Seitz et al. 2001). As the densities of these prey decrease, the proportional mortality decreases, so the prey have a low-density refuge from predation. The functional response can be changed both qualitatively and quantitatively by changes in habitat (Lipcius and Hines 1986, Dittel et al. 1995, Long et al. 2012b, predator size and density (Mansour andLipcius 1991, Gonza´lez-Sua´rez et al. 2011), the presence and density of alternative prey (Chesson 1989, Long et al. 2012a, and the spatial arrangement of prey (Hines et al. 2009, Long andHines 2012). Blue crabs have a type III FR toward M. balthica (Eggleston et al. 1992, Taylor andEggleston 2000), but as hypoxia should decrease predator search time for M. balthica, we hypothesized that it may change the FR from a type III to a type I or II. This could have negative implications for the population, because both type I and type II FRs can lead to the local extinction of prey populations.

Survey of populations in the York and Rappahannock Rivers
Macoma balthica populations in the York and Rappahannock Rivers were surveyed by box-coring before and after summer hypoxic events in 2003 and 2004. We randomly selected 10 sites within each of three depth strata: shallow (,3 m), mid-depth (3-6 m), and deep (.6 m). In 2003 in the Rappahannock, five sites were sampled in each stratum. We stratified by depth to ensure that we sampled enough hypoxic sites; hypoxia is strongly correlated with depth in both rivers (Kuo and Neilson 1987 We used a 25 3 25 cm Gray O'Hara box core to take benthic samples on each date. In 2003, from the center of each box core, one 10 cm diameter core was taken down to 15 cm and sieved on a 0.5-mm screen. The remainder of the core was sieved on a 2-mm screen. Everything retained on the sieves was frozen and processed in the laboratory. Bottom temperature, salinity, DO, water depth, and box-core penetration were measured. A small number of cores had shallow penetration due to hard substrate; cores with penetration ,15 cm were not used for 2-mm analyses, and cores with penetrations ,5 cm were not used for 0.5-mm analyses Comtois 1985, Long and. In 2004, the sampling was the same, except that we took three replicate box cores at each site and only sub-cored the first of these. The 2-mm and 0.5-mm samples were sorted in the laboratory, and all M. balthica were removed and measured. The 0.5-mm samples were stained with Rose Bengal before sorting and were used to quantify newly recruited M. balthica (shell length 3-5 mm). Juveniles (shell length 5-18 mm) and adults (shell length . 18 mm) were retained and enumerated from the 2-mm samples. In 2004, the dry and ash-free dry masses (AFDM) were determined for all M. balthica, except for the recruits.
We used maximum likelihood fitting to examine the effect of DO on the survival of M. balthica over the summer in 2004. We calculated survival (S ) at each site as the average density after hypoxia divided by the average density before hypoxia. We modeled survival using the following equation: where S is survival (dependent variable), and DO is average summer DO (independent variable). The parameters are S max (maximum survival), DO 50 (DO concentration at which survival is 50%), and b (curvature). We assumed a normal distribution of errors with a standard deviation (SD) that increased with DO: SD ¼ sDO 3 , where s is a constant. We examined this and a second model that allowed survival to change with the factor river (MATLAB, Version 7.2.0.232; MathWorks 2006). Akaike's information criterion, corrected for small sample size (AIC c ), which is a measure of the relative parsimony of a model, was used to select the best (i.e., most parsimonious) out of a set of models (Burnham and Anderson 2002). In this and subsequent analyses, we considered models with DAIC c , 2 to be plausible (Burnham and Anderson 2002). M. balthica shell length frequencies from 2004 were analyzed using maximum likelihood techniques (MAT-LAB, Version 7.2.0.232; MathWorks 2006). We hy-pothesized a priori that the shell lengths could be modeled as a bimodal normal distribution, with two means and two standard deviations (representing the adult and the juvenile cohorts) and a mixing parameter. Juvenile M. balthica recruit in the fall and spring in the Chesapeake (Shaw 1965). Female M. balthica start producing eggs at a shell length of ;18 mm (Long and Seitz 2008b), and only some will attain that size by the fall of their first year (see data in Results: Survey populations in the York and Rappahannock Rivers); most will not reproduce until their second year when they join the adult size class.
We considered three factors: river, time (before or after hypoxia), and DO. As we were interested in sublethal effects of DO, we used a DO cutoff above hypoxia, dividing the sites into low-DO sites (average summer DO , 3.5 mg/L) and high-DO sites (average summer DO . 3.5 mg/L). This cutoff provided an adequate sample size for the low-DO sites and represents the point where survival is reduced by 15% of the maximum as estimated by our field survival analysis (in Methods: Survey populations in the York and Rappahannock Rivers). Thus, clams at sites with DO . 3.5 mg/L are unlikely to be stressed by low DO, whereas clams at sites with DO , 3.5 mg/L are expected to be stressed (Long et al. 2008b). We considered nine models where the mean size of each cohort and the mixing parameter were linear functions of these factors (Table 1). Time was considered important because mean size was expected to increase over the summer as the clams grew. In addition, time was expected to affect the mixing parameter because the smallest clams should have higher mortality rates than larger clams. We also hypothesized that low DO would decrease average size of clams because stressed individuals have less energy to devote to growth. We also considered that sites with low DO would have populations of clams that were dominated by juveniles, as mortality from low DO the previous year would reduce the adult segment of the population. Thus, we allowed DO levels to influence the mixing parameter. The potential effects of river were unknown, but we considered that each river would have a different effect on mean size due to differing environmental conditions. We calculated AIC c values for each model and used them to select the best-fitting model. We then calculated 95% confidence intervals using the slice method for each of the parameters for the best-fitting model.
We used maximum likelihood techniques to determine effects of time, river, and DO on the length-mass relationship in M. balthica (MATLAB, Version 7.2.0.232; MathWorks 2006). We used a standard exponential growth curve as the base model: We assumed the errors had a normal distribution, with the standard deviation increasing with the Length 3 . We analyzed 10 a priori models wherein b was a linear function of time, river, and DO (Table 3). Again, time was considered important because clams were expected to lose mass in the fall after spawning. Low DO was expected to decrease mass, and river was included to account for differing environmental conditions. We also considered an interaction term between time and DO, as the effects of low DO were likely only to be significant after summer hypoxia. The AIC c values were calculated as before (see Methods: Survey populations in the York and Rappahannock Rivers) for each model, and the 95% confidence intervals calculated by the slice method for the parameters in the best-fitting model.

Effects of hypoxia on growth and reproduction
We performed a laboratory experiment to test the effects of sublethal hypoxia on the growth and reproduction of M. balthica. We used 16 40-L plastic tanks, which were filled with 20 cm of mud collected from shallow coves in the York River where M. balthica are abundant, and that was passed through a 2-mm sieve to remove bivalves. M. balthica with shell lengths .10 mm were collected from muddy sites in the York River and brought back to the laboratory. Only clams with whole shells and a quick siphon-withdrawal reflex were used. Twenty clams were randomly selected for each based on Akaike's information criterion with small sample size correction (AIC c ).
tank. Each clam's shell was dried, measured, and marked with a letter using permanent ink. Clams were allowed to bury for at least 24 h before the experiment was started to allow acclimation to normal burial depth (Seitz et al. 2003b). Any clams that failed to bury in 24 h were removed and replaced.
We produced hypoxic water by flowing unfiltered seawater through a fluidized mud-bed reactor (Long et al. 2008b). Hypoxic water was stored in a 400-L distribution tank, which was bubbled with nitrogen to maintain hypoxia. A second distribution tank was filled with unfiltered seawater that was kept normoxic by vigorously bubbling air through it. Each tank was randomly assigned to one of four treatments: (1) hypoxic (DO , 2.0 mg/L), (2) moderately hypoxic (DO 2.5 6 0.5 mg/L), (3) nearly normoxic (DO 3.5 6 0.5 mg/L), and (4) normoxic (DO . 4.0 mg/L). To enable the measurement of DO, temperature, and salinity without opening the tank and letting oxygen in, each tank was covered with a lid that had a hole drilled in it; a cork was inserted into the hole when measurements were not being taken.
The DO in tanks was reduced to the nominal DO range by dispensing differing proportions of hypoxic and normoxic waters from the distribution tanks into the experimental tanks. Temperature, DO, and salinity were monitored approximately twice a day, and the flow of water to each tank was adjusted when necessary to keep each tank within its nominal DO range. The experiment ran for 28 d, from 21 August to 18 September 2006, after which all the tanks were returned to normoxic conditions. Tanks were monitored and a high flow rate of normoxic water was maintained until 10 November 2006, when the clams were harvested by sieving the mud through a 5-mm screen. Mortality rates were calculated for each tank using the equation S ¼ Ne Àmt where S is the number of survivors, N the initial number, t the time in days, and m the mortality rate. The final shell length was measured and the growth of each clam determined as the change in shell length. All clams were opened, gender was determined (male, female, or juvenile), and the AFDM of each clam was measured.
For fecundity measurements, all female clams were opened and homogenized in 15 mL of phosphatebuffered saline (PBS; 10 mmol/L Na 2 HPO 4 , 150 mmol/ L NaCl, pH 7.2), and then centrifuged at 16 671.3 m/s 2 (1700 g) for 1 h. An enzyme-linked immunosorbent assay (ELISA) using a monoclonal antibody specific to an HSP70 protein that is indicative of oogenesis in M. balthica (Bromage et al. 2008) was performed on each sample according to the protocol in Long et al. (2008a) to quantify egg production. Clams from a hypoxic tank and a normoxic tank were used to establish the relationship between the concentration of mb-HSP70 and egg concentration (Long et al. 2008a). Samples of eggs were removed from the ovaries, suspended in PBS, and counted with a hemocytometer. The samples were centrifuged at 147 099.8 m/s 2 (15 000 g) for 1 h at 48C, homogenized and centrifuged again. The mb-HSP70 concentration in the supernatant was determined by ELISA, and the protein concentration was determined with a bicinchoninic acid protein assay kit (Sigma) using bovine serum albumin for the standards. We calculated the amount of protein and mb-HSP70 in each egg. The mb-HSP70 : egg ratio was used to calculate the number of eggs produced by each clam.
Temperature, growth, and fecundity were analyzed using a nested analysis of variance (ANOVA) with tank nested within DO treatment. The protein : egg and mb-HSP70 : egg ratios were analyzed using a nested ANOVA with clam number nested within DO treatment. Mortality rates were analyzed with a one-way ANOVA using DO treatment as the factor.

Hypoxia, density dependence, and predation
In summer 2006, we performed a caging experiment to determine the effects of hypoxia upon predation on M. balthica. We randomly selected four sites at each of two depths: shallow (3-4 m) and deep (10-12 m) in the York River near Gloucester Point (Fig. 1). Sites were selected to keep environmental factors, especially sediment type, temperature, and salinity, consistent across the sites. The experiment was performed twice, once during 6-29 June 2006, and again during 3-31 August 2006. During the June experiment, both the deep and shallow sites remained normoxic. During August, the deep sites experienced episodic hypoxia and the shallow sites remained normoxic. At each site, four 50 3 50 cm plots were established by SCUBA and marked with a frame made of polyvinyl chloride (PVC) pipe. Each plot was a combination of a caging treatment (caged, partially caged [to control for caging artifacts], or uncaged) and a clam density treatment (high [40 clams/plot ¼ 160/m 2 ] or low [10 clams/plot ¼ 40/m 2 ]). These densities are within the normal range in the York River (Seitz et al. 2005). We used these two densities because it is possible to distinguish between type II and type III FRs using only two prey densities, one low and the other moderately high, when the FR has been previously determined (Seitz et al. 2001); a type II FR will lead to higher predation rate at low densities, whereas a type III FR will lead to the opposite (Taylor and Eggleston 2000). The plots were as follows: (1) caged, high density; (2) partially caged, high density; (3) uncaged, high density; (4) uncaged, low density. Live M. balthica of shell length .10 mm were collected from the York River and brought back to the laboratory, where they were marked with permanent ink and allowed to recover for 24 h. Only healthy clams with whole un-chipped shells, and a quick siphon-withdrawal response were used. They were transplanted in each plot at the nominal density and covered with a full cage for at least 24 h to allow undisturbed burial and acclimation (Seitz et al. 2001). The full cage was then removed from the uncaged plots, replaced with a partial cage in the partially caged plots, and left in place in the caged plots. The cages were made from galvanized steel hardwire cloth (1-cm mesh), and were 50 3 50 cm 3 14 cm high. Partial cages had a 25 3 25 cm square cut in the center of the top, and a 25 cm wide 3 7 cm tall hole cut in the center of each side to allow access by predators. All cages were inserted 7 cm into the sediment so that the bottoms of the holes in the sides of the partial cages were flush with the sediment surface.
Each experiment lasted ;28 d before all plots were resampled with a suction sampler to a depth of 40 cm (Eggleston et al. 1992). Samples were brought back to the laboratory and frozen before processing. Each sample was sorted, marked M. balthica were counted, and the percentage of recovery was calculated for each plot. Marked, whole M. balthica shells were counted as recovered when we calculated predation because they represented non-predatory mortality . Unmarked bivalves were identified and measured.
A continuous water quality recorder (model Hydrolab DS5X with a luminescent DO probe; Hach Environmental, Loveland, Colorado, USA) was deployed at one of the deep sites during both experiments. This instrument recorded temperature, salinity, and DO every 5 min and was downloaded and serviced every 7-10 d. During the first deployment in June, the probe temporarily failed. A five-point running average was applied to the data to smooth it. Spot measurements were made throughout the experiment at each of the sites using a DO probe (YSI Model 85; Yellow Springs Instruments, Dayton, Ohio, USA).
Predators were sampled by the Virginia Institute of Marine Science's (Gloucester Point, Virginia, USA) Juvenile Finfish and Blue Crab Trawl Survey. The survey samples sites throughout the Virginia portion of Chesapeake Bay and its tributaries. Monthly, 5-min trawls are performed using a 9-m semi-balloon otter trawl with a 38.1-mm stretch-mesh body and a 6.35-mm mesh cod-end liner. All animals caught in the trawl were identified to species and length was measured.
Predation rate was calculated using N ¼ N 0 e Àpt , where N is the number of clams surviving after time t, N 0 is the initial number, and p is the rate of predation . The rate of predation may include a caging effect (from the partial cages) and a density-dependent effect. We calculated the predation rate for each of the uncaged high-density, uncaged lowdensity, and partially caged plots using the number of clams recovered in the caged plot at that site as N to account for recovery efficiency (i.e., non-recovered clams). For the low-density plots, we assumed that recovery efficiency was not density dependent, and used the proportional recovery 3 10 as N. In some instances, we had cage failure; either the cage was broken or missing when we returned to the site, most likely caused by fisheries activity. In these cases, we used the average recovery in the caged plots within that time and depth as N. We also estimated non-predatory mortality in the uncaged and partially caged plots with the same equation using the number of dead, whole M. balthica recovered as ''1 À S,'' and the number recovered in the caged plot as N. We analyzed predation rates using an ANOVA with depth (deep or shallow), time (June or August), and plot (uncaged high density, uncaged low density, or partially caged) as main factors, and site (nested within depth and time) as a blocking factor. Where an interaction term was significant, a Student-Newman-Keuls post hoc multiple comparison test (SNK test) was performed.
Modeling the M. balthica population Density-independent model.-We used the results from the previous experiments to parameterize and analyze a matrix model for the M. balthica population in the York River (Caswell 2001). We used two age classes, adults (A, reproductively mature, .1 year old) and juveniles (J, ,1 year old), and divided the metapopulation into two populations: Individuals in areas that remain normoxic (subscript N), and those in areas that go hypoxic (subscript H), giving the following population vector: The population census took place after recruitment in the spring, before predators entered the river, when the abundance is greatest. We assume that predation only occurred from the late spring through early fall when predators are common in the system and actively foraging (Hines et al. 1990), which we estimated as 1 May to 30 September (153 d). For the density-independent model, we assumed that predation was density independent, as the results of our predation experiment suggested. Survival over the summer in normoxic areas was calculated as: S N ¼ e Àð pNþmNÞT ; where S is the proportion surviving, p is the predation rate, m is the non-predatory mortality rate, and T is the time in days from the start of the summer. To calculate survival in hypoxic areas, we assumed that, during normoxia, predatory and non-predatory mortality are the same in all areas of the river. Thus, survival in hypoxic areas is given by S H ¼ e Àð pHþmHÞTH e Àð pNþmNÞðTÀTHÞ , where T H is the time in days that the river experiences episodic periodic hypoxia, which we estimated as 90 d (Kuo and Neilson 1987). We estimated p from our predation experiment and m from our laboratory experiments for both hypoxic and normoxic areas. Our initial estimate of p N resulted in summertime survival rates that were much lower than those observed, so we lowered p N within the 95% CI until the predicted S N was similar to that observed in the field. Our estimate of p N may have been high because our predation experiments were performed during periods of high predator abundance and did not account for lower predation rates in the early spring or late fall. The predicted S H was similar to observed survival, so we did not adjust the values of p H . Fecundity was estimated from our laboratory experiment for clams in both hypoxic and normoxic areas. We assumed that recruitment R was equal to the number of eggs produced times a constant survival term. We calculated R as the average density of recruits in the spring of 2004 divided by average density of adult clams in fall 2003. Juvenile clams that attain a shell length of 18 mm will reproduce in their first year and produce about 20% as many eggs as an adult (Long et al. 2008a). We estimated the proportion of juveniles that attain this shell length (M ) in hypoxic and normoxic areas by the predicted shell length frequencies calculated from the population survey. We assumed that recruitment was equal across all areas, and that d is the proportion of the river that remains normoxic and (1 À d) is the proportion that goes hypoxic. We further assumed that overwinter survival is 100%, that all juveniles graduate to the adult size class after their first year, and that clams do not move from one area to another after recruitment. From this, we constructed the population projection matrix (Eq. 1).
This model contains three parameters that allowed us to explicitly examine effects of a changing hypoxic regime on the population. The duration of hypoxia, t H , and the areal extent of hypoxia, ''1 À d,'' are explicit in the model. We calculated the dominant eigenvalue and eigenvectors, and performed sensitivity and elasticity analyses on the matrix. We further examined the effects of various parameters on the long-term stability of the population and on the population's stable age structure.
Incorporating density dependence.-The predator functional response (FR) of blue crab Callinectes sapidus to M. balthica in Chesapeake Bay is type III, or sigmoid, in the laboratory (Eggleston et al. 1992) and is the same type as the guild FR in the field (Seitz et al. 2001). We used this to create a hybrid model that contains both a continuous part (mortality and growth over the summer) and a discrete part (reproduction and graduation of size classes). We used the same population vector as above (see Methods: Density independent model ), changing the units to density (clams/m 2 ). This gave rise to the population projection matrix (Eq. 2), similar to Eq. 1, but lacking survival terms.
The following differential equation describes summer mortality: where N is the density of M. balthica, m is non-predatory mortality, g(N ) is the predator FR, and P is the density of predators. We assumed a type III FR (Hassell et al. 1977): where b and c are parameters of the instantaneous attack rate (Hassell et al. 1977), and T hand is the handling time. We used parameters estimated by Eggleston et al. (1992) for C. sapidus feeding on M. balthica in mud: b ¼ 0.006, c ¼ 0.001, and T hand ¼ 4.06 h. Predator densities in the York River are strongly correlated with prey densities (Seitz et al. 2003a. We used the relationship between blue crab density and clam density in the York River described by Seitz et al. (2003a) to estimate the predator density at each clam density: The rate of predation during hypoxia was increased by a factor of 3.0, the increase in the predation rate observed in our predation experiment. We calculated the number and biomass of clams consumed by predators and lost to nonpredatory mortality by integrating these ordinary differential equations: where N P and N m are the number of clams/m 2 that suffer predatory and non-predatory mortality, respectively; B P and B m are the biomasses of those clams; and B c is the biomass of an individual clam and is a function of the length, L. We used the length-mass relationship determined above (see Methods: Survey populations in the York and Rappahannock Rivers) for clams in the York River in the spring of 2004 for B c . L(t) is the von Bertalanffy (1938) growth function: L(t) ¼ L ' À (L ' À L 0 )e Àkt , where L ' is the maximum size, L 0 is the initial size, and k is the growth rate. We estimated these parameters using preliminary results from a markrecapture experiment in the York River (L ' ¼ 32.7 mm; k ¼ 0.85 yr À1 ; sample size ¼ 441; B. Brylawski, unpublished data). L 0 was the average length of clams in each age class and area from the York River in the spring of 2004. The census was taken immediately after recruitment in the spring. We then ran the continuous portion for five months (153 days), after which the population vector was multiplied by the population project matrix. We assumed that hypoxia took place in the middle of the summer. We performed sensitivity analysis on the model and examined how changing the hypoxic duration and extent affected the equilibrium populations of M. balthica, and how it affected the number and biomass of clams consumed by predators and lost to nonpredatory mortality. We further examined the intensity of hypoxia (i.e., the duration of hypoxic episodes and the proximity of DO levels to anoxia) with the m H parameter, which should increase with increasing hypoxic intensity.

Survey of populations in the York and Rappahannock Rivers
Hypoxia was more extensive in the Rappahannock River than in the York River, and more extensive in 2004 than 2003 (Fig. 1a The Macoma balthica populations in both rivers followed a general pattern where recruitment pulses in the fall and spring (in Chesapeake M. balthica recruits twice a year; Shaw 1965) were followed by decline in density over the summer (Fig. 2). Densities diminished in all areas over the summer, but in hypoxic areas they declined to local extinction (Fig. 2). Recruitment was patchy, but no clear trends differentiated recruitment in hypoxic and normoxic sites. In February 2004, recruitment tended, though not significantly, to be higher in hypoxic sites than in normoxic sites in the Rappahannock (Fig. 2); in May 2004, normoxic sites in the York had higher densities than hypoxic sites (Fig. 2).
Survival was increased with average summer DO (Fig.  3). Average survival was given by the following model: Including the river parameter increased the AIC c by 2.5, indicating that survival did not differ between the rivers (Burnham and Anderson 2002).
In the best-fit model for shell length frequencies, the means of both distributions were functions of time, river, and DO, and the mixing parameter was a function of time and DO (Table 1). Under low DO, mean shell length of juveniles (l 1 ) increased, whereas it decreased in adults (l 2 ), and the population was shifted toward juveniles by 20% (Table 2, Fig. 4). Juveniles grew by ;14 mm over the course of the summer (i.e., between the before-hypoxia and the after-hypoxia samples) and the adults remained essentially unchanged, while the populations shifted toward the adults by 40% (Table 2, Fig. 4). River had a mixed effect: In the York, the average shell length of the juveniles was slightly larger, but that of adults was slightly smaller (Table 2, Fig. 4).
Two models of the relationship between shell length and AFDM had substantial weight (Table 3). Both of these indicated an effect of time and river; clams in the fall (after hypoxia) had, on average, lower AFDM at a given shell length than in spring (before hypoxia), and clams in the York had a higher AFDM than those in the Rappahannock (Fig. 5). The evidence for an effect of low DO (Tables 3 and 4) was weak; the parameter represented a biologically insignificant change of ;0.25% in b. Thus, the results presented in Fig. 5 do not include this term; while the length-mass relationship changed seasonally, it did not vary between areas that experienced high and low DO.

Effects of hypoxia on growth, mortality, and reproduction
The DO for each experimental tank in hypoxia treatments generally remained within the nominal range (Fig. 6), but it varied over time, which is characteristic of DO in natural systems such as the York River. Average temperature during the experiment was 24.28 6 1.18C (mean 6 SD), which did not differ among hypoxia treatments (ANOVA, F 3, 687 ¼ 1.36, P ¼ 0.25) or tanks (ANOVA, F 12, 687 ¼ 0.67, P ¼ 0.79).
Clam growth was minimal over the experimental period and averaged 0.02 6 0.38 mm shell length. Mortality did not differ among treatments during the experiment (ANOVA, F 3,12 ¼ 1.12, P ¼ 0.38); average mortality rate was 0.0014 6 0.00025 clams/day, which was equivalent to an average of 10.6% mortality over the duration of the experiment.
Hypoxia decreased the number of eggs produced by female M. balthica, but increased the protein investment in each. The ratio of mb-HSP70 : eggs was influenced by hypoxia treatment (ANOVA, F 1,12 ¼ 7.28, P ¼ 0.019) and clam number (ANOVA, F 3,12 ¼ 11.33, P ¼ 0.001), with the hypoxic tanks having a higher ratio than normoxic tanks ( Fig. 7a; Tukey's test). The same pattern was observed in the ratio of protein : egg; both hypoxia treatment (ANOVA, F 1,12 ¼ 9.40, P ¼ 0.01) and clam (ANOVA, F 3,12 ¼ 7.45, P ¼ 0.004) had an effect, with the hypoxic tanks having a higher ratio than normoxic tanks (Tukey's test; Fig. 7b). As the mb-HSP70 : egg ratio differed between the hypoxic and normoxic treatments, we were able to calculate egg production for clams from hypoxic and normoxic treatments, but not from the two intermediate treatments. We therefore performed two analyses. We ran a nested ANOVA on the concentrations of mb-HSP70 in each clam (a measure of maternal investment), with Treatment and Tank nested within Treatment as factors, using all female clams from all treatments. We then calculated the number of eggs produced by each clam in the normoxic and hypoxic tanks and performed a second identical nested ANOVA. Neither Treatment (ANOVA, F 3,56 ¼ 0.31, P ¼ 0.818; Fig. 7c) nor Tank (ANOVA, F 11,56 ¼ 1.05, P ¼ 0.417) affected the concentration of mb-HSP70. Treatment (ANOVA, F 1,31 ¼ 4.72, P ¼ 0.038), but not Tank (ANOVA, F 6,31 ¼ 1.31, P ¼ 0.28) affected egg production. Clams in hypoxic tanks produced 40% fewer eggs than those in normoxic tanks (Tukey's test; Fig. 7d).

Hypoxia, density dependence, and predation in the field
During the June field experiment, DO measurements remained well above 2.0 mg/L in both the deep and shallow areas, as expected (Fig. 8a). Although we lost some data from the water quality recorder at the beginning of the experiment, we were able to collect data during neap tides, when hypoxia is most likely (Kuo and Neilson 1987). In August, the shallow sites remained normoxic during the experiment, but the deep sites suffered from one strong hypoxic episode that lasted about five days, beginning 3 August, and also had several short excursions into hypoxia starting 28 August (Fig. 8b). During both experiments, temperatures in the deep area were ;18C cooler than shallow areas and ;1 psu more saline.
Hypoxia increased the rate of predation, but not of non-predatory mortality. Average recovery of clams in caged plots was 92% in June and 85% in August. In June, we lost the caged, uncaged high, and partially caged plots from one deep site, and the uncaged low plot from a second deep site in August. This left us with at least three replicates of each plot type in each depth and time period, which was sufficient for analysis. There was a significant interaction effect between time and depth (ANOVA, F 1,19 ¼ 14.9, P ¼ 0.001). Predation rate in the deep, hypoxic sites during August was significantly higher than all other depths and time periods (SNK, P , 0.01; Fig. 9a). There was also a significant effect of site on the predation rate (ANOVA, F 11,19 ¼ 2.49, P ¼ 0.02); however, the only difference was between the two deep sites in August that had lost one or more plots, and therefore, had sample sizes of 1 and 2, making the comparison unreliable (Tukey's test). There was a significant interaction effect of time and depth on nonpredatory mortality (ANOVA, F 1,19 ¼ 4.7, P ¼ 0.043). The rate of non-predatory mortality was higher in the deep sites in August than in the deep sites in June (SNK, P , 0.05), but there was no difference between the deep and shallow sites in August.
Density did not have a significant effect on predation rate, but the trends are interesting (Fig. 9b). One of the major difficulties in the interpretation of the results was the presence of ambient unmarked clams that increased the actual density of M. balthica above the desired density and blurred the distinction between nominally high-and low-density plots. This problem was especially prominent at shallow sites. In the deep sites during hypoxia, predation rates were so high that we recovered no live clams in one of our partially caged, high-density plots. With such high losses in all density plots, differences in proportional mortality were slight. In August, however, during hypoxia, proportional mortalities in the high-density plots at shallow sites tended to be higher, whereas proportional mortalities in lowdensity plots at deep sites tended to be higher.
The overall abundance of predators was low in June-July, and was dominated by hogchoker and Atlantic croaker (Fig. 10). In August-September, the abundance of predators increased 10-fold and the relative number of blue crabs increased, though the population was still dominated by piscine predators (Fig. 10).

Model simulation
Density-independent model.-Our population projection matrix, based on the parameter estimates (Table 5) The dominant eigenvalue (k) of this matrix is 1.04, indicating a growing population (k . 1 implies that the population is increasing, and k , 1 that population is decreasing; Caswell 2001). Sensitivity analysis indicated that the dominant eigenvalue was sensitive to reproduction and survival of clams in normoxic areas, as well as recruits produced by maturing juveniles in hypoxic areas that recruit into normoxic areas (Fig. 11a). Elasticity analysis showed that the three matrix elements that contribute substantially to k are survival of clams in normoxic areas, and recruitment into normoxic areas from adults in normoxic areas (sum of these elements ¼ 0.98). Clams in hypoxic areas contributed little to k (sum of elements including hypoxic clams ¼ 0.006; Fig. 11b).
The dominant eigenvalue, k, was sensitive to changes in the extent of normoxia, d, along its full range and was relatively insensitive to changes in the duration of hypoxia, T H , between about 153 and 40 days, but increased rapidly between 40 and 0 days (Fig. 12). As predicted from the elasticity analysis, k was almost completely insensitive to reproduction from clams in hypoxic areas (Fig. 13); k was fairly sensitive to changes in predation in hypoxic areas, p H, from 0 to 0.15 d À1 and insensitive at larger values (Fig. 13).
The stable age structure associated with the densityindependent matrix showed that ;80% of the population was juveniles and ;20% was adults in normoxic areas; adults in hypoxic areas were a negligible part of the population (Fig. 14). The stable age structure was not much affected by changes in t H or p H . At small values of t H and p H , A H approached A N ; this is expected, because when conditions in hypoxic areas approach those in normoxic areas, the age structure in hypoxic areas should resemble that in normoxic areas (Fig.  14a, d). As d decreased, the population became progressively dominated by J H , and as R N decreased, the population became dominated by A N (Fig. 14c, d).
Density-dependent model.-The model predicted cyclical population dynamics with recruitment in the spring and fall followed by mortality over the summer (Fig. 15), which is typical of M. balthica (Holland et al. 1987). The population reached equilibrium rapidly, after about 10-25 generations. In normoxic areas in the early spring, adult density increased to ;200 clams/m 2 , when the previous year's juveniles became adults, while juvenile density was near 400 clams/m 2 , which is a moderate density of clams in shallow areas in the York . Summer mortality was much higher in hypoxic areas, and the total population dropped to ,12 clams/ m 2 , compared with 190 clams/m 2 in normoxic areas (Fig.  15). Numerical analysis suggested that, under certain conditions, there are two stable equilibria, one high and one low density, and an unstable equilibrium at zero. Which equilibrium the population approaches depended on the initial population vector. The population also FIG. 5. Observed (points) and predicted (lines) relationship between shell length and ash-free dry mass (AFDM) for M. balthica in the (a) York and (b) Rappahannock Rivers in 2004. Notes: The base model is Mass ¼ aLength b , which assumes a normal distribution of errors with a standard deviation ¼ sLength 3 . Parentheses indicate that the parameter is a linear function of the parenthetical factors. See Table 1 for clarification of abbreviations.
switched from the high-to low-density equilibrium if population densities dropped below a critical density (which represents a second unstable equilibrium), as might be caused by some catastrophic event. Under other conditions, the lower equilibrium became the only stable one, regardless of the initial population.
The model predicts that the metapopulation of M. balthica is sensitive to the spatial extent (1 À d ), and duration (t H ) of hypoxia, but not to its intensity (m H ; Fig. 16). The density of clams decreased with increasing extent of hypoxia, exhibiting a threshold response at d ¼ 0.49 (when just over half the river is hypoxic), below which the metapopulation was functionally extinct (Fig.  16a). Additionally, as the extent of hypoxia increased, the threshold density (i.e., the density at which the population will switch from the high-to the lowequilibrium densities) increased, indicating that the resilience of the population to perturbation decreased.
The population of adult clams in hypoxic areas decreased with increasing hypoxic duration until it became functionally extinct when the duration was greater than 100 days (Fig. 16b). Varying the intensity of hypoxia did not change the equilibrium densities of the population even when the rate of non-predatory mortality was increased to that observed under anoxia (Fig. 16c).
Hypoxia also affected predation and trophic transfer (Fig. 17). Increasing the extent of hypoxia decreased the number and the biomass of clams eaten in both hypoxic and normoxic areas (Fig. 17a, b). Although the number of clams consumed in hypoxic areas was generally predicted to be higher than that in normoxic areas, the biomass was generally lower because a greater proportion of those clams were juveniles. The exception for this trend occurred when the duration of hypoxia was ,40 days, in which case, the model predicted that the   Tables 1 and 3  biomass consumed in hypoxic areas would be greater than that in normoxic areas (Fig. 17c). The amount of biomass lost to non-predatory mortality was low and almost always much lower than that consumed by predators. However, increasing the intensity of hypoxia, that is, increasing the rate of non-predatory mortality in hypoxic areas, changed this ratio. When the rate of predation and the rate of non-predatory mortality were approximately equal, the biomass lost to each was equal, but as the non-predatory mortality rates approached those observed under anoxic conditions, the biomass lost to non-predatory mortality greatly exceeded that lost to predation (Fig. 17d).

Effects of hypoxia at the individual level
Low DO has a multitude of effects on individual M. balthica. Periodic hypoxia, such as was experienced in our lowest DO treatment in our tank experiment and in August in the deep area during our field predation experiment, did not increase the non-predatory mortality rate, as seen in other studies (Modig and Ó lafsson 1998). Periodic hypoxia may also induce a physiological tolerance to longer term hypoxia (Altieri 2006). Only after long periods of continuous hypoxia did nonpredatory mortality in M. balthica increase in laboratory experiments (e.g., de Zwaan et al. 2001, Seitz et al. 2003b). In our studies, sublethal hypoxia decreased egg production in M. balthica, similar to other species (e.g., Chrysaora quinquecerrha, Condon et al. 2001; Mnemiopsis leidyi and Chrysaora quinquecirrha, Grove and Breitburg 2005; Acartia tonsa, Richmond et al. 2006). Interestingly, in our experiment, as the number of eggs decreased, their protein content increased. This is consistent with predictions based on evolutionary and life-history theory that individuals in stressful environments will maximize their fitness by producing fewer offspring, but investing more in each (McGingley et al. 1987, Sibly et al. 1988. M. balthica also have behavioral responses to hypoxia that include increasing siphon extension (Seitz et al. 2003b) and reducing burial depth (Brafield 1963, Long et al. 2008b) to gain access to higher DO levels above the benthic boundary layer. Both of these responses take place within 24-48 hours of the initiation of hypoxia.
Our analysis of shell length frequencies supports the hypothesis that hypoxia affects the growth of M. balthica, as is common for other species (e.g., McNatt and Rice 2004, Grove and Breitburg 2005, Richmond et al. 2006. Hypoxia can also result in delayed reproduction (Richmond et al. 2006, Brouwer et al. 2007). In May, the juveniles in areas that experienced low DO had shell lengths that were ;1 mm longer than those in high-DO areas; however, the adult populations in the same area had shell lengths that were ;1 mm shorter. This suggests that the negative effects of low DO on growth are subtle and only detectable after a year, explaining why we did not discern an effect in our laboratory experiment; the duration was not long enough for differences in growth to be noticeable. Juveniles may be larger in hypoxic sites because of lower ambient clam densities in these areas (from reduced abundances due to hypoxic episodes the previous year) precluding intraspecific competition for food (Rhodes andYoung 1970, Woodin 1974); M. balthica growth, though not mortality, is affected by density (Ó lafsson 1986). Low DO did not have a substantial effect on the relationship between clam length and mass, indicating that the clams use energy for the maintenance of optimal mass rather than for growth. The difference in the length-mass relationship between the samples taken before and after hypoxia are due to seasonal differences possibly related to reproductive status.

Effects of hypoxia at the population level
Combined, the multiple sublethal effects of hypoxic stress on M. balthica can lead to temporary, local extinction of populations, a substantial decrease in production, transfer of biomass up the food web, and possibly the collapse of the metapopulation. In the York and Rappahannock Rivers, this combination of effects leads to source-sink dynamics along a DO gradient. In our field study, mortality was complete in hypoxic areas, leading to local extinction. Thus, areas where DO drops below ;1.5 mg/L represent a black-hole sink (sensu Gomulkiewicz et al. 1999) for M. balthica, and areas where the DO averages between ;1.5 and 2.5 mg/L represent a traditional sink habitat. In a sink habitat, the rate of reproduction is not sufficient to maintain the population in that area (Pulliam 1988, Dias 1996, Holt 1997. In a black-hole sink the organisms recruit, but die before they can reproduce (Dias 1996, Gomulkiewicz et al. 1999). Thus, M. balthica that recruit into hypoxic areas will not be able to maintain the population, as they produce fewer eggs, or will die before reproduction. Though source-sink dynamics likely occurs in various marine invertebrate species (e.g., Lipcius et al. 1997, Lipcius and Ralph 2011, ours is the first empirical evidence for source-sink dynamics in an infaunal bivalve. The populations of M. balthica in the York and Rappahannock Rivers do not appear to suffer from recruitment limitation; recruitment is strong both in areas that remain normoxic as well as in areas that experience hypoxia. Indeed, recruitment in the Rappahannock River in the spring of 2004 tended to be higher in the areas that typically go hypoxic, possibly because of reduced competition (Rhodes andYoung 1970, Woodin 1974) due to death of the macrobenthos during the previous year's hypoxia. This pattern of strong recruitment of M. balthica in areas that go hypoxic occurs in Chesapeake Bay (Holland et al. 1987) and elsewhere (Powers et al. 2005). This combination of high recruitment and low survival in hypoxic areas gives rise to a population dominated by juveniles, which is both what we observed, and what was predicted by both of our models. In areas between 1.5 and 2.5 mg/L DO, low survival of M. balthica, combined with reduced growth and fecundity, result in a level of reproduction insufficient for maintenance of the population (e.g., Kreuzer and Huntly 2003). The populations in these areas are replenished each year by recruits from shallow, normoxic habitats that act as sources.
Our models of the population dynamics of M. balthica predict what may happen to the population under changing hypoxic conditions. For example, source-sink population dynamics can lead to destabilization of the population as a whole. Our model predicts that when the extent and duration of hypoxia reach a threshold, the population of M. balthica will decline toward extinction, as more general models predict (e.g., Delibes et al. 2001b). Increasing the intensity of hypoxia, or decreasing reproduction in clams in hypoxic areas, has little effect on population stability, because mortality in FIG. 10. Total abundance of benthic predators in the York River during the summer of 2006. Points to the left of the vertical line are plotted against the left y-axis and those to the right against the right y-axis. A change of scale was necessary to show the details of the predator guild during the middle of the summer. Arrows indicate the time periods during which the first and second predation experiments (Exp.) were performed. Predators are: summer flounder, Paralichthys dentatus; spot, Leiostomus xanthurus; hogchoker, Trinectes maculatus; croaker, Micropogonias undulatus; and blue crab, Callinectes sapidus. Notes: Parameters are: p N and p H are the predation rates in normoxic and hypoxic areas, respectively; m N and m H are the non-predatory mortality rates in normoxic and hypoxic areas; T is the duration of summer; T H is the duration of hypoxia; M N and M H are the proportions of the spring recruits that reach reproductive maturity in normoxic and hypoxic areas; R N and R H are the number of recruits produced per mature female in normoxic and hypoxic areas; and d is the areal proportion of the river that remains normoxic throughout the summer. hypoxic areas is already near 100%. This is why the dominant eigenvalue of the model is sensitive only to changes in survival and reproduction in normoxic areas. Our model assumes that the system is closed, and that M. balthica larvae do not immigrate from areas outside the coupled source-sink pair; such immigration could help stabilize a population. However, if the population in a given river required external inputs of larvae from other rivers to maintain stability, then the population in that river would be functioning as a sink. Increasing the scale in this way does not affect the general prediction that, once threshold conditions are reached on the largest effective scale, the meta-population will decline.
For species that have a more active role in choosing their habitat, a sink habitat is often occupied by individuals forced out of source habitats by competition (Pulliam 1988). In this case, sinks are only populated after the source habitats have reached carrying capacity, and thus, the total population size with the sink habitat is greater than it would be with the source habitat alone (Dias 1996). However, for species that disperse passively, or are unable to distinguish between source and sink habitats, a sink habitat can reduce the population in the source (Pulliam and Danielson 1991, Delibes et al.  . R H is expressed as a proportion of the reproduction of clams in normoxic areas so that when R H ¼ 1, they reproduce at the same rate. Note that because mortality is density independent, the effect of varying p H is identical to varying m H , the non-predatory mortality rate. 2001a, b). This latter case is what our density-dependent model predicts, and what our observations support. Increasing the extent and duration of hypoxia in the model can cause the equilibrium population density in the normoxic areas to drop by .40%, and the population density of M. balthica in the Rappahannock, which generally has a greater extent and duration of hypoxia, is lower than in the York, as our model predicts. This suggests that hypoxia may contribute to the lower M. balthica population densities in the Rappahannock.
Our model further predicts that there are alternative stable states for the M. balthica population, and that increasing hypoxic stress can decrease the resilience of the population by increasing the threshold density at which a switch between the high-and low-density states occurs. A switch in states is predicted to cause the population to decline by two orders of magnitude or more. A decline in the density of a species of this magnitude has already occurred for Crassostrea virginica, the eastern oyster, in both Chesapeake Bay and the Neuse River, where declining water quality combined with high fishing mortality, habitat degradation, and disease have led to population collapse (Rothschild et al. 1994, Lenihan andPeterson 1998).

Effects of hypoxia on trophic interactions
Our predation study demonstrated that hypoxia can have a substantial effect on trophic transfer. The rate of predation in the deep areas in August, during hypoxia, was increased threefold compared with predation rates under normoxic conditions. During the experimental period, there was a substantial hypoxic episode at the beginning of the experiment, and a slight one at the end. The duration of the first episode was five days, which is long enough for M. balthica to exhibit a behavioral response, but not long enough for non-predatory mortality to occur (Seitz et al. 2003b, Long et al. 2008a). These results corroborate those obtained in a similar experiment the previous year , providing strong evidence that hypoxic episodes of sublethal duration considerably increase vulnerability of M. balthica to predation, and that predators in the system take advantage of this opportunity. A similar behaviorally induced change in prey behavior in response to episodic hypoxia can increase predation rates on bivalves in other systems Mistri 2011, 2012).
Although motile predators typically avoid hypoxic areas (Pihl et al. 1991, Eby and Crowder 2004, there are examples of predators entering hypoxic water to forage both in estuarine (Nestlerode andDiaz 1998, Roberts et al. 2001) and freshwater systems (Rahel and Nutzman 1994), indicating that predation can occur during hypoxic events. In the Gulf of Mexico, cownose rays (Rhinoptera bonasus) will feed in extensive hypoxic areas, but spend most of their time swimming above hypoxic areas in shallower, normoxic waters (Craig et al. 2010). In addition, predators do not avoid episodic hypoxia to the same extent as they do chronic hypoxia (Bell et al. 2003, Bell and. Finally, M. balthica will exhibit a behavioral response to hypoxia at DO concentrations higher than the concentrations that piscine and crustacean predators avoid , Long et al. 2008b, Vaquer-Sunyer and Duarte 2008. Epibenthic predators were present in abundance in the general area of our field experiments and may have taken advantage of stressed infaunal prey during or just after hypoxic events. Taken together, these responses offer a large window of opportunity for increased predation during and after hypoxia. In previous studies in the York River, predation rate was lower during hypoxia (Nestlerode and Diaz 1998), so the most likely mechanism for the FIG. 16. Effects of (a) the proportional extent of hypoxia (1 À d ), (b) the duration of hypoxia (T H ), and (c) the intensity of hypoxia as measured by the rate of non-predatory mortality (m H ) on the equilibrium density of M. balthica in the York River, predicted by the density-dependent model. Each point represents the predicted density after 50 generations immediately after spring recruitment. Note that, because we assume recruits settle equally in all areas, the density of juveniles (J ) in normoxic (subscript N) and hypoxic (subscript H) areas is the same, whereas the density of adults (A) varies between areas. Effect of the extent of hypoxia (1 À d ) on the threshold density (the density of A N at which the population will switch between the high and low equilibria, given a starting density of 170 clams/m 2 in all other ages and areas) is shown in plot (a).
increased predation in hypoxic areas in our study is reinvasion by predators immediately after a hypoxic event (Pihl et al. 1991), before prey species have had a chance to recover (Pihl et al. 1992). If this is the case, then short-lived hypoxic episodes with a duration of ;1 week are important for trophic transfer.
Although we did not detect a significant effect of density on predation rates, our results suggest that hypoxia changes the functional response of predators to M. balthica. In shallow areas, ambient clams increased the density in our plots above our nominal density, making detection of a density-dependent response difficult. However, the trend in shallow plots of higher rates of predation in higher density plots is consistent with a type III functional response. This response has been demonstrated for M. balthica in multiple field and laboratory experiments (Mansour and Lipcius 1991, Eggleston et al. 1992, Seitz et al. 2001, as well as for other bivalves (Lipcius and Hines 1986). In our deep sites, the predation rate was so high that we observed complete local extinction in one of our high-density plots, which is why high-and low-density plots did not differ significantly, though the trend was for higher predation rates in low-density plots. This trend is consistent with a type II functional response and in opposition to a type III functional response (Taylor and Eggleston 2000). Finally, a type III functional response implies the existence of a low-density refuge at which prey are safe from predation (Eggleston et al. 1992). For M. balthica, this refuge exists between 12 and 50 clams/ m 2 in the field in the York River (Seitz et al. 2001). Our observation of extinction in a high-density plot, and the fact that our low-density hypoxic plots were reduced from an average of 40 clam/m 2 to 8 clams/m 2 , shows that the clams no longer have a low-density refuge from predation under hypoxia, providing further evidence against type III and support for a type II functional response.

Effects of hypoxia at the ecosystem level
Whereas the distinction between predatory and nonpredatory mortality is not of great importance at the individual and population levels, at an ecosystem level it is critical (Altieri andWitman 2006, Altieri 2008). Clams that die from non-predatory mortality enter the microbial loop and are lost to higher trophic levels (Baird et al. 2004), but clams that are consumed by predators are passed up the food chain (Pihl et al. 1992). If the majority of the biomass that is lost during hypoxia goes to predators, then the overall effect on the ecosystem may be beneficial in terms of trophic transfer up the food web, in contrast to the situation when the biomass enters the microbial loop directly. In the latter situation, breakdown of the biomass by microbes likely aggravates hypoxia when microbes respire during metabolism (Altieri and Witman 2006).
In the York River, predation rate increased under hypoxic conditions, but the non-predatory mortality did not, as compared to normoxic sites. This increases both the rate of trophic transfer and the ratio of biomass consumed to biomass lost to microbes; predators consumed so many clams that there were few left to die from other causes. Our model predicts that many more clams are consumed in hypoxic areas than in normoxic ones. However, because the population of clams in hypoxic areas is dominated by juveniles, the biomass consumed is generally substantially less in hypoxic areas.
The magnitude of hypoxia is important in determining the rate of trophic transfer. As the areal extent of hypoxia increases, population levels in both the hypoxic and normoxic areas decrease. This lowers the number of clams available and subsequently the number and biomass of the clams consumed by predators in all areas. The duration of hypoxia has an interesting range of effects on trophic transfer. When the duration is low (,40 days), biomass consumed in hypoxic areas is higher than in normoxic areas. A short duration of hypoxia affords predators a narrow window of opportunity for higher consumption rates, accounting for the high amount of biomass consumed during these periods. However, the level of predation is not high enough to cause local extinction of the population, so benthic production remains high, and some of the juveniles survive to become adults, which increases the available biomass. Indeed, short-term hypoxia can actually create a partial refuge from predation if, as is the case for the hard clam Mercenaria mercenaria, (1) the hypoxia is not severe enough to increase prey mortality, (2) the hypoxia is severe enough to reduce predator abundance, and (3) the hypoxia does not cause the prey to exhibit a response that increases venerability to predation (Altieri 2008). In our model, as the duration of hypoxia increases, predators consume greater proportions of the prey, production is negatively affected, and the total biomass consumed decreases. The situation is analogous to human overharvesting of natural resources leading to a decrease in total harvest; maximum yields occur at a moderate rate of exploitation (e.g., Rieman and Beamesderfer 1990). Similarly, intermediate disturbance leads to highest density and diversity of benthic species (Connell 1978, Sommer 1995. Increasing the intensity of hypoxia (i.e., moving from episodic to chronic hypoxia and decreasing the DO concentrations towards anoxia) magnifies the rate of non-predatory mortality. This substantially decreases the biomass of clams consumed by predators and increases the biomass lost to the microbial loop. At mortality rates equal to those we observed for M. balthica under anoxia, most biomass in hypoxic sites enters the microbial loop, as predicted by our and other models (Baird et al. 2004), and as observed in systems with more intense hypoxia (Powers et al. 2005, Altieri andWitman 2006).
An important ecosystem-level effect of hypoxia is the decrease in filtration by the benthic assemblage. When M. balthica densities decrease in association with hypoxia, the overall density and biomass of the benthic assemblage, which includes deposit, filter, and suspension feeders, also declines. Consequently, filtration rate will decline along with the biomass, because filtration rate is proportional to the biomass of suspension feeders. M. balthica is a facultative suspension and deposit feeder (Kammermans 1994). Filtration by bivalves and other benthic fauna has the potential to reduce the effects of eutrophication in the system by removing labile organic material from the system (Officer et al. 1982, Newell et al. 2007, so reduction in bivalve biomass by hypoxia could, through a feedback mechanism, exacerbate the problems caused by eutrophication (Altieri and Witman 2006), as has been proposed for oysters in Chesapeake Bay (Newell et al. 2007). Although the presence of clams can increase nutrient concentrations in water (Doering et al. 1987), either through excretion (Nakamura et al. 1988) or bioturbation (itself another important ecosystem service provided by clams that can increase oxygenation of sediments and benthic-pelagic coupling; Jones et al. 1994), thus potentially decreasing or negating the positive effects of filtration (Chan and Bendell 2013), the net effect is likely to be removal of nutrients (Nakamura et al. 1988) and an increase in local water quality (e.g., Wall et al. 2011).
Another potential ecosystem-level effect of hypoxia is an increased density of predators in normoxic areas due to flight from hypoxic waters (Eby and Crowder 2002). This can result in an increase in predation in shallow habitats (Lenihan et al. 2001), and potentially increased injuries and cannibalism within predator populations , Aumann et al. 2006). In addition, aggregations of organisms along the edge of hypoxic patches may increase their vulnerability to predation or fishing pressure (Craig andCrowder 2005, Craig 2012). We did not observe an increase in predation in the shallow areas during hypoxic episodes, suggesting that habitat compression was not occurring in this system (Pihl et al. 1992). This may be either because the areal extent of hypoxia was not great enough to produce an observable effect, because the predator or prey density had changed between the experiments, or because predators were inefficient at avoiding episodic hypoxia (Bell et al. 2003, Bell and. These ecosystem-level effects are sobering, especially if one considers that the effects of hypoxia on benthic species are the same as for pelagic species, where hypoxia similarly affects growth (Breitburg 1992, Grove andBreitburg 2005), reproduction (Condon et al. 2001, McNatt and Rice 2004, Richmond et al. 2006, mortality (Breitburg 1992), and predation (Breitburg et al. 1994(Breitburg et al. , 1997. Adding the effects of hypoxia on benthos and the water column further decreases the overall services the ecosystem is able to provide, and decreases its stability and resilience (Baird et al. 2004).

CONCLUSIONS
In this study, we have explored the effects of hypoxia on the biomass-dominant species in an estuarine system and, through modeling, integrated those effects over time, space, and the ecosystem. The effects cascade up through the various scales. Individual effects on growth, reproduction, and mortality decrease benthic abundance, biomass, and production in local populations. Further, behavioral changes in individual M. balthica result in a higher rate of predation, increasing trophic transfer, but decreasing benthic biomass and secondary production. Integration over space, by considering linkages between subpopulations, gives rise to sourcesink population dynamics and results in the population in hypoxic areas reducing the equilibrium abundances in normoxic areas. We predict that a short duration of hypoxia may result in a slight increase in ecosystem services (e.g., increased trophic transfer and production in commercially important species), but over the long term, increasing the extent and duration of hypoxia could lead to a population crash and elimination of these important ecosystem services. Out of these dynamics emerge ecosystem-level effects, where shortterm increases in predation may lead to long-term decreases in energy transfer from primary producers to predators, and a feedback mechanism between the lost benthic filtration and the development of hypoxia could further decrease the stability and resilience of the system to eutrophication. These results strongly suggest that efforts to control or reduce eutrophication and control or reduce hypoxia are warranted (Rabalais et al. 2007).