Central to ecology and resource management is knowledge of the spatiotemporal scales at which demographic rates vary and the ecological consequences of demographic variation, such as that due to density dependence. We quantified the spatiotemporal variation in eastern oyster Crassostrea virginica recruitment, density, growth, and survival and assessed density dependence within a network of no-take reserves in Pamlico Sound, North Carolina. From 2006 to 2008, average oyster recruitment and total density increased 15- and fivefold, respectively. The unprecedentedly high oyster densities in certain reserves (up to 6,500/m2 at the end of the study) modified demographic rates such that further density increases were regulated by density-dependent survival. Oyster demographic rates varied significantly among reserves at relatively small spatial scales (20 km). Certain reserves were the strong “recruiters,” others the fast “growers,” and yet others the high “survivors.” Cohort dynamics altered the demographic rank order such that the demographically “best” reserves varied intra- and interannually. From a management perspective, the prevalence of density-dependent survival suggests that the oysters in this system are habitat rather than recruitment limited, which may minimize the utility of stock enhancement programs. Addition of habitat (i.e., artificial reefs) should focus on reserves characterized by high recruitment but density-dependent growth and survival. This study (1) supports the efficacy of marine reserves in rapidly increasing the density and age—size structure of protected species, (2) highlights the need for spatially explicit demographic data to support multifaceted management objectives, and (3) when combined with evidence of reserve larval connectivity, provides support for applying metapopulation concepts to this reserve system.
The dynamic nature of marine populations is driven by spatiotemporal variation in several demographic rates, including recruitment, growth, survival, and reproduction (Gotelli 2001; Jennings et al. 2001). Many marine populations are demographically open, whereby recruitment, is uncoupled from local reproduction (Caley et al. 1996). At sufficiently large spatial sales, local populations connected by larval dispersal form a demographically closed metapopulation (Levins 1969; Hanski and Gilpin 1991). Inherent in the metapopulation concept is the incorporation of data on spatiotemporal variation in local demographic rates that can contribute disproportionately to metapopulation dynamics (Figueira 2009). Knowledge of metapopulation dynamics is critical to the management of marine systems (Sale et al. 2006), particularly with the growing use of marine protected areas. In this work, we quantify spatiotemporal variation in local demographic rates of eastern oyster Crassostrea virginica within a network (i.e., metapopulation; Haase et al. 2012) of no-take broodstock reserves.
Population density, a key response variable in assessing a species' abundance patterns and certain types of resource management strategies, can modify demographic rates and regulate populations (Caley et al. 1996). For population regulation to occur, at least one demographic rate (e.g., growth or survival) must be density dependent such that the demographic rate decreases as population density increases (Murdoch 1994). In the absence of density dependence, populations are considered recruitment Subject editor: Anne Hollowed, Alaska Fisheries Science Center, Seattle limited (sensu Doherty 1982; see reviews by Caley et al. 1996; Hixon et al. 2002). The potential for (meta)population size to be regulated by density dependence is a central tenet in ecology, with important implications for managing fisheries and their habitats (Sánchez Lizaso et al. 2000; Jennings et al. 2001; Gust 2004; Chockley et al. 2008; Figueira 2009; Lorenzen et al. 2010).
Marine protected areas and, more specifically, no-take marine reserves where all forms of extractive activities are prohibited, are a potentially powerful restoration, management, and conservation tool (Plan Development Team 1990; Roberts and Polunin 1993; Allison et al. 1998; Halpern and Warner 2002; Lubchenco et al. 2003; Eggleston and Parsons 2008; Gaines et al. 2010). Demographic benefits within reserve boundaries are well documented (Halpern and Warner 2002; Halpern 2003; Claudet et al. 2008; Lester et al. 2009; Babcock et al. 2010). In a global meta-analysis of 124 marine reserves using three functional taxonomic groups, Lester et al. (2009) reported that relative to controls outside the reserve or before reserve establishment, organism density and body size (among other biological responses) within reserves were, on average, 166% and 28% higher, respectively. Perhaps just as noteworthy was the variability in biological response to reserve designation; changes in reserve density ranged from approximately -90% to 2000% (mean = 166% increase; median = ∼50% increase; Lester et al. 2009).
Clearly, demographic variability is to be expected across global scales and taxonomic groups—the scales of metaanalyses. What about at the smaller spatial scales (10–100 km) over which marine reserve networks and metapopulations often operate? Are demographic patterns consistent over space and time such that certain reserves are demographically the “best” on a consistent basis? Or do reserve networks consist of a shifting mosaic of reserve demographic rank order that can serve as a buffer to biotic and abiotic variability (Gaines et al. 2010)? In this study, we address these questions as they relate to the success of a large-scale eastern oyster restoration effort consisting of a network of marine reserves.
The global decline of native oyster populations due primarily to long-term overfishing and habitat destruction (Gross and Smyth 1946; Kirby 2004; Beck et al. 2011) has negatively impacted estuarine ecosystems (Jackson et al. 2001). Worldwide, restoration efforts are under way to restore populations of these economically (Rothschild et al. 1994; Mackenzie 2007) and ecologically (Wells 1961; Ulanowicz and Tuttle 1992; Coen et al. 1999; Peterson et al. 2003; Rodney and Paynter 2006; Fulford et al. 2010) important bivalves. Restoration techniques range from constructing 3-dimensional artificial reefs within reserves to hatchery-based stock enhancement (Coen and Luckenbach 2000; Laing et al. 2006; Quan et al. 2009; Beck et al. 2011). Recent efforts to restore eastern oysters (referred to simply as oysters hereafter) along the U.S. Atlantic coast have focused on establishing and assessing the efficacy of reserves (Powers et al. 2009; Schulte et al. 2009; Paynter et al. 2010). Fundamental to the success of such restoration efforts are ecological processes such as recruitment, growth, survival, and density dependence (Lenihan et al. 1999; Coen and Luckenbach 2000; Mann and Evans 2004; Brumbaugh and Coen 2009; Powers et al. 2009; Paynter et al. 2010).
In this article we present the results of a 3-year study that quantified the spatiotemporal variation in demographic rates of six oyster reserves in Pamlico Sound, North Carolina. Specifically, we (1) quantified oyster density, recruitment, growth, and survival at each reserve, (2) compared demographic rates over space and time, and (3) investigated the potential for density dependence in space (with demographic rates related to differences in density among reserves within a cohort) and time (with demographic rates related to differences in density among cohorts within a reserve).
STUDY SITE
The Croatan—Albemarle—Pamlico estuarine system (CAPES) is the largest lagoonal system and second largest estuarine system in the United States, covering an area of ∼6,600 km2 (Epperly and Ross 1986; Luettich et al. 2002; Figure 1). The CAPES is separated from the Atlantic Ocean by the Outer Banks barrier island chain, which limits exchange with shelf waters to four narrow inlets (Pietrafesa et al. 1986; Lin et al. 2007). Pamlico Sound, the largest component of the CAPES (∼120×40 km; Figure 1), provides ideal nursery habitat for many estuarine-dependent finfish and shellfish populations (Epperly and Ross 1986; Eggleston et al. 2010) due to its shallow water depths (mean, ∼4.5 m), high primary productivity, relatively stable salinities, and low currents and tidal amplitudes (Paerl et al. 2001).
FIGURE 1.
Map of the Croatan—Albemarle—Pamlico estuarine system, with the locations of oyster reserves in Pamlico Sound indicated by squares (squares not to scale). The labeled reserves—West Bay (WB), Ocracoke (OC), Hatteras (HA), Crab Hole (CH), Bluff Point (BP), and Deep Bay (DB)—were studied.
TABLE 1.
Oyster reserves in Pamlico Sound sampled from 2006 through 2008. The number of mounds refers to those created in 2003 (i.e., those sampled with quadrats) unless otherwise noted. In 2006 and 2007, the number of quadrat replicates was allocated to be proportional to the number of mounds created during 2003. In 2008, three quadrat replicates were sampled at all reserves. The area of hard substrate within each reserve was estimated from side-scan sonar surveys (Ballance and Eggleston 2008). Salinity (psu) was measured during sampling events. Mean salinity (psu) was averaged among years within months. For reserve locations, see Figure 1.
Since 1996, the North Carolina Division of Marine Fisheries has established 10 subtidal oyster reserves in meso- and polyhaline (salinity, 13–29 practical salinity units [psu]) waters of Pamlico Sound where oyster harvest and the use of bottomdisturbing fishing gear is prohibited (Table 1; Figure 1). Within reserve boundaries, high-relief (∼2 m off the bottom), coneshaped mounds were constructed from ∼150 tons of ∼0.5-m limestone riprap to provide larval settlement substrate. Recent research in this reserve system indicated that oyster densities at 4 of the 10 reserves ranged from 0 to 97/m2 in 2002–2003 (Powers et al. 2009) and that potential larval connectivity between reserves was present (15 of the possible 90 inter-reserve connections) yet asymmetrical, with some reserves providing more connections than others (Haase et al. 2012). In this work, we chose 6 of the 10 oyster reserves for study—West Bay, Ocracoke, Hatteras, Crab Hole, Bluff Point, and Deep Bay (Figure 1)—to span the length—width axis of Pamlico Sound in an attempt to capture the spatial variation in oyster reserve demographic rates. The reserves studied ranged in area from 0.03 to 0.19 km2, the distance between them from 20 to 105 km, and age from 1 to 10 years (Table 1).
METHODS
Demographics: density and length frequency.—Oyster density was estimated at each of the six oyster reserves (Figure 1) during May—June and July—August 2006–2008 using scubabased 0.25-m2 quadrat sampling and hand excavation from randomly selected (in June 2006) reef mounds within each reserve. In 2006 and 2007, quadrat sampling effort was proportional to the number of reefs constructed during 2003, when present (Table 1), to minimize the confounding effects of reef age on de mographic rates. In 2008, three quadrat replicates were surveyed at all reserves. Oysters, riprap, and shell material were excavated to a depth of ∼15 cm within each quadrat (sensu Powers et al. 2009), brought to the surface in mesh bags, and processed within 24 h of collection. All live oysters were counted. The left valve length (LVL), i.e., the distance from the umbo region of the shell to the anterior shell margin, was measured with calipers to the nearest 0.1 mm on a subsample (≥ 1/8) of the excavated oysters. Length frequencies and counts from each quadrat were scaled to oysters/m2.
Replicate quadrat length frequency distributions were compared using a two-sample Kolmogorov—Smirnov test with Bonferroni adjustments for multiple comparisons. Only 5 of the 124 quadrat replicates were significantly different (P ≤ 0.02; otherwise P ≥ 0.08) from their counterparts within the same reserve on the same sampling date, so the quadrat replicates at each reserve for a given sampling period were combined. A modal analysis of oyster length frequency distributions (replicates combined) was conducted to follow the progression of cohort LVL and density at each reserve. Size-class intervals were specified at 5 mm. Modes in each monthly length frequency distribution were identified using NORMSEP (FiSAT), which treats length frequency distributions as a mixture of normal distributions and applies a maximum likelihood procedure to separate the length frequency distributions into their normal components. Each mode was characterized by a mean LVL, SD, and density. The number of distributions (i.e., modes) was specified a priori at one and increased in number until the separation index of modal means was <2 (∼2 SD; Jennings et al. 2001 and references therein). Modes were assumed to represent distinct cohorts, and cohorts were assumed to progress to the next larger mode in subsequent length frequency distributions (Table 2; Appendix 1). Based on chi-square goodness-of-fit tests, modal decomposition provided adequate fits to the observed length frequency data in 25 of 30 instances (P ≥ 0.06; Table 2; Appendix 1).
TABLE 2.
Modal mean ± SD length and density (N/m2) estimates from length frequency analysis of three oyster cohorts—June 2006 (J06), August 2006 (A06), and May 2007 (M07)—at West Bay (WB), Ocracoke (OC), Hatteras (HA), Crab Hole (CH), Bluff Point (BP), and Deep Bay (DB) oyster reserves in Pamlico Sound. Chi-square goodness-of-fit tests are reported for all cohorts (not only those shown here). Lack of fit is indicated by significant P-values (shown in bold italics). When cohorts could not be decomposed, a single modal mean, SD, and N are reported. See Appendix 1 for length frequency data and fitted modal means.
We tested whether mean oyster density (total, recruit [≤25 mm], and legal size [>75 mm]) varied significantly among six sampling dates (June 2006, August 2006, May 2007, August 2007, May 2008, and July 2008) and six reserves (West Bay, Ocracoke, Hatteras, Crab Hole, Bluff Point, and Deep Bay), with a two-way analysis of variance (ANOVA). For all ANOVAs, each quadrat sample was considered a replicate. In cases in which there was a significant sampling date × reserve interaction, we used a one-way ANOVA to test for differences in mean oyster density among reserves within a sampling date and among sampling dates within a reserve (i.e., simple effects). Post hoc pairwise comparisons were performed on all significant simple effects, and Bonferroni adjustments were used to determine the significance of pairwise contrasts.
Demographics: growth and survival.—To quantify the spatiotemporal variation in oyster growth and survival in the absence of the confounding effects of local variation in oyster density and reef architecture, we employed five replicate standardized oyster settlement substrates (trays) at each reserve from June 2006 to October 2008 as part of a field mark—recapture study. Each open-topped plastic mesh tray (∼0.25 m2) was uniquely numbered and standardized with 40 disarticulated oyster shells cemented in a 1-in layer of concrete. To allow settlement of natural oyster cohorts, scuba divers deployed trays in June 2006, August 2006, and May 2007 near the top of five randomly selected (in June 2006) reef mounds at each reserve prior to primary and secondary soundwide settlement peaks, which typically occur during June and August, respectively (Eggleston et al. 2011). Settlement peaks can also occur in September and October, particularly along the Outer Banks (Eggleston et al. 2011), but we did not monitor these cohorts. We followed the demographic fate of the three oyster cohorts that initially colonized the settlement trays, i.e., the cohorts from June 2006, August 2006, and May 2007. Cohort age was determined by the date of settlement (on which age = 0 d), which was assumed to occur the day of tray deployment (giving us conservative growth estimates). The trays were retrieved by scuba divers at bimonthly intervals from May to October during 2006–2008. Upon retrieval, the trays were digitally photographed alongside a metric ruler for georectification. Prior to being photographed, ∼20 oysters per tray were marked with nail polish for identification in the photographs and measured with calipers to calibrate the measurements of LVL from image analysis (see below). Fouling organisms were removed from the trays to ensure a clear photograph (sensu Roegner and Mann 1995; Bishop and Peterson 2006). Afterwards, divers returned the trays to their original locations.
Image analysis software (Image Pro-Plus) was used to estimate oyster growth and survival from the time series of digital photographs. The initial photograph in each time series (i.e., cohort) for each settlement tray was used to enumerate all individuals for density estimates. Forty individuals in each cohort (when present) were haphazardly chosen, “marked,” measured, and subsequently tracked through time to quantify their growth and survival. Survivorship was quantified as the proportion of the original cohort that survived to the start of each age/census in the time series (Gotelli 2001). Individuals were assumed to have survived at each census if they were present in photographs with both valves and minimal valve gape. If they were present and alive during a census, the LVL of tracked individuals in each cohort was then measured to the nearest 0.1 mm. Due to the three dimensional (i.e., vertical) nature of oyster growth and two dimensional digital images, length corrections were applied to oysters growing at angles between 30° and 90°. Two separate linear regressions were used to correct (i.e., inverse regression) the image analysis estimates of LVL for oysters oriented from 30° to 75° and 75° to 90°. Following inverse regression, predicted LVL was highly correlated (r2 > 0.92) with caliper LVL along the 1:1 line (i.e., there was no bias; slope and y-intercept of residuals = 0: t < 0.01, P > 0.99).
To model the size and survival at age for each cohort at each reserve, nonlinear least-squares regression was used to fit cohort- and reserve-specific seasonalized (i.e., winter growth stasis) von Bertalanffy functions (VBF) and Weibull functions (Pinder et al. 1978), respectively. These two models were selected to describe oyster growth and survivorship due to their flexibility and ecologically meaningful parameters. Left valve length (mm) at age t, LVL(t), was calculated according to the equation
where LVL∞ is the maximum asymptotic left valve length, K is a curvature parameter that determines how fast LVL∞ is approached, t0 is the theoretical age (years) when LVL is zero, C is related to the magnitude of seasonal oscillations, and tw is the time at which growth is slowest. The parameters t0 and C were constrained to values ≥ 0.0 and 0.7–1.0, respectively. The probability of survival to age t, S(t), was calculated according to the equation
where b is the scale parameter, and c is the shape parameter. Estimates of c > 1, = 1, and <1 imply type I (the mortality rate increases with age), type II (the mortality rate is constant), and type III (the mortality rate decreases with age) survivorship curves, respectively.
To test for spatial and temporal variation in growth and survivorship, nested comparisons of seasonalized VBF and Weibull functions were conducted among cohorts within a reserve and among reserves within a cohort using likelihood ratio tests (Kimura 1980). We first tested for differences in modeled growth and survivorship in time among all cohorts within a reserve across all model parameters. To maximize sample sizes, tests among cohorts were conducted over the complete age range of each cohort as opposed to a common age range. When growth and survivorship models were significantly different in these global tests, we (1) tested whether the model differences among cohorts were caused by any single parameter, (2) used the results of the tests in (1) to select the most parsimonious models, and (3) conducted pairwise comparisons between cohorts within a reserve using the models selected in (2). We applied the growth and survivorship models selected in (2) to conduct global tests in space among reserves within a cohort across all model parameters. When the global test was significant, we repeated steps 1–3 for comparisons among reserves. In single-parameter tests of the growth models, t0 and C were not included because these parameters were largely constrained (see above). Bonferroni adjustments were used to determine the significance of multiple pairwise comparisons.
Density dependence.—We tested for density-dependent oyster growth and survival in space (data pooled among reserves within a cohort) and time (data pooled among cohorts within a reserve) on standardized (settlement trays) and nonstandardized (quadrat) substrates. Analysis of covariance (ANCOVA) models were used to test for density dependence. Model terms included LVL at age 1 or the proportion surviving to age 1 as the response variables, cohort or reserve as the categorical factors for spatial and temporal tests, respectively, age-0 recruit density or total density (for settlement trays, recruit density = total density) as covariates, and the interaction between covariate and categorical factors (sensu Steele and Forrester 2005). The interaction term was removed from models when nonsignificant; density dependence was determined by a significant covariate in the reduced model. When the interaction term was significant, separate slope models were fitted to the relationship between the response variable and recruit or total density; slopes significantly different from 0 were interpreted as indicating density dependence (sensu Connell 1985).
Age-0 recruit density and total density were estimated from settlement trays and quadrats as explained above in the sections titled “Demographics: density and length frequency” and “Demographics: growth and survival.” Age-1 survival and mean LVL were estimated differently for trays and quadrats. For settlement trays, the proportion surviving to age 1 and mean LVL at age 1 were estimated from a subset of oysters within each tray that were censured as explained in “Demographics: growth and survival.” For quadrats, modal analysis, as explained in “Demographics: density and length frequency,” provided mean cohort LVL and cohort density over time. Survival was calculated as the proportion of the recruiting cohort remaining at each census. We assumed that August recruits settled in May or June to estimate LVL and survival at age 1.
RESULTS
Density
A total of 55,868 oysters were counted and 20,862 measured for LVL from the six oyster reserves over the six sampling periods from June 2006 to July 2008. Mean total, recruit, and legal oyster density varied significantly by reserve (two-way ANOVA: F5, 96 ≥ 17.1, P < 0.0001) and time (F5, 96 ≥ 19.1, P < 0.0001); however, a significant reserve × time interaction (F25, 96 ≥ 4.4, P < 0.0001; Figure 2) precluded contrasts across the main effects.
Total oyster density.—Mean total oyster density pooled among reserves increased by 451%, from 686 to 3,782 oysters/m2 during the 3-year study period (Figures 2a, 3; Appendix 2). Within reserves, the percent increase in total oyster densities over time ranged from 87% (Hatteras) to 1,662% (Bluff Point (Figure 3). During this time, the age and size structure of oysters within reserves changed from predominately bimodal to quadrimodal (Appendix 1). Mean total oyster density increased significantly over time at all reserves (F5, 118 ≥ 15.9, P < 0.0001) except Hatteras and Deep Bay (F5, 118 ≤ 1.9, P ≥ 0.1). Significant increases in total oyster density occurred, on average, at annual intervals concurrent with August peaks in recruitment (Figure 2a; Appendices 1, 2). Comparisons of mean total oyster densities among reserves were significant at each point of time (F5, 118 ≥ 2.4, P ≤ 0.04) except August 2006 (F5, 118 = 1.3, P = 0.3). The Deep Bay reserve consistently ranked lowest in total oyster density at each point of time, whereas the Crab Hole and Ocracoke reserves typically ranked highest (Appendix 2).
Recruit density.—Oyster recruit (≤25 mm LVL) densities estimated from settlement trays and quadrats were correlated along the 1:1 line (r = 0.87; slope and y-intercept of residuals = 0; t ≤ 2.1, P ≥ 0.06). Mean oyster recruit density pooled across reserves increased by 1,480% during the 3-year study period (Appendix 2). Oyster recruitment occurred at all reserves and in all sampling periods with the exception of Deep Bay in June 2006 (Figure 2b; Appendices 1, 2), when mean recruit density was homogenous among all reserves (F5, 118 = 1.0, P = 0.5). Oyster recruitment was ∼3 times greater during July—August (May—June cohort) than during May—June (September—October cohort). Mean oyster recruit density varied significantly over time at each reserve (F5, 118 ≥ 5.1, P ≤ 0.0004) with the exception of Hatteras and Deep Bay (F5, 118 ≤ 1.6, P ≥ 0.2). Mean oyster recruit density groupings were nearly identical to total oyster densities beginning in May 2007, with the Deep Bay reserve consistently ranking the lowest and the Crab Hole reserve the highest at each time (Appendix 2). Legal density.—The mean density of legal-sized (>75 mm LVL) oysters increased, on average, from 40 to 300/m2 and by ≥240% at each reserve over time (Figure 2c; Appendix 2). The density of legal-sized oysters at West Bay, Ocracoke, Hatteras, and Bluff Point increased over time (F5, 118 ≥ 3.0, P < 0.01). Beginning in May 2007, the mean density of legal oysters varied among reserves at each time (F5, 118 ≥ 9.9, P < 0.0001). Ocracoke consistently ranked the highest in density of mean legal-sized oysters. Unlike its position with respect to total and recruit oyster densities, Crab Hole consistently ranked lowest in density of legal oysters (Figure 2c; Appendix 2).
FIGURE 2.
Mean (a) total, (b) recruit (≤25 mm LVL), and (c) legal (>75 mm LVL) oyster density at the six reserves studied (see Figure 1) at various times. Note the differences in the scale of the z-axis. See Appendix 2 for density values.
Growth and Survival
A total of 2,706 oysters from the June 2006, August 2006, and May 2007 cohorts were tracked on settlement trays until October 2008 for growth and survival analyses (Table 3). To test for a settlement tray artifact, we compared tray- and quadrat-based estimates of mean cohort LVL and survivorship. Estimates of LVL and survivorship were correlated between sampling methods (LVL: r = 0.96, P < 0.0001; survivorship: r = 0.52, P < 0.001; Appendices 3, 4), although tray-based estimates of LVL and survivorship were, on average, 9% and 18% higher than the corresponding estimates from quadrats. The differences in estimated demographic rates between the two methods were significant for LVL at Crab Hole (t = 5.5, df = 14, P < 0.0001) and for survivorship at West Bay and Deep Bay (t ≥ 3.2, df ≥ 9, P ≤ 0.005). In general, there was little evidence for a settlement tray artifact whereby trays biased growth and survivorship estimates.
Growth.—Seasonalized VBF accurately described the size-at-age relationship for each oyster cohort at all reserves (R2 ≥ 0.84; Appendix 3). Age-specific oyster growth, as modeled by seasonalized VBF, was significantly different in time among cohorts within a reserve (χ2 ≥ 19.5, df = 5, P < 0.002). Temporal differences in modeled growth among cohorts were driven by K; LVL∞ and tw were invariant among cohorts. The minimum and maximum model-predicted size-at-age among cohorts within a reserve varied by < 1–57%, 7–40%, and 2–30% at 0.5, 1.0, and 2.5 years of age, respectively. Cohorts from all reserves excluding Crab Hole and the May 2007 cohort at Ocracoke and Hatteras reached legal size—76 mm—by age 1.5 years (Figure 4; Appendix 3). No reserve exhibited significant differences in modeled growth among all three cohorts.
The seasonalized VBF, in which LVL∞ and tw were set equal among cohorts within a reserve, differed spatially among the reserves within a cohort (Χ;2 ≥ 44.2, df = 5, P < 0.0001; Figure 4). Spatial differences in modeled growth were driven by differences in LVL∞ and K (χ2 ≤ 10.9, df = 1, P ≤ 0.001) but not tw (χ2 ≤ 3.0, df = 1, P ≥ 0.08). Across reserves, estimates of LVL∞ for the June 2006, August 2006, and May 2007 cohorts fell in the ranges 67.7–107.4, 81.2–185.3, and 61.2–87.3, respectively (Appendix 3). Likewise, estimates of K were highly variable across reserves within each cohort—June 2006: 0.9–1.9, August 2006: 0.4–1.1, and May 2007: 1.4–2.2 (Appendix 3). Among reserves, the minimum and maximum model-predicted size at age within a cohort varied by 24–42%, 19–31%, and 31–36% at 0.5, 1.0, and 2.5 years of age, respectively (Figure 4). The cohorts at Crab Hole were consistently the smallest at a given age (Figure 4). The patterns were less consistent among other reserves where size-at-age rank order varied among cohorts. For example, Ocracoke ranked in the largest grouping for the June 2006 and August 2006 cohorts but in the smallest grouping for the May 2007 cohort (Figure 4).
TABLE 3.
Initial information for three oyster cohorts tracked through time for growth and survival on settlement trays at six oyster reserves in Pamlico Sound. The settlement date is assumed to be the date of tray deployment.
Survival.—Age-specific oyster survivorship was accurately modeled by the Weibull function for each reserve—cohort combination (R2 ≥ 0.70) with the exception of the May 2007 cohort at Deep Bay, where model fits were poor (R2 = 0.44; Appendix 4). Age-specific oyster survivorship, as modeled by Weibull functions, was significantly different in time among cohorts within a reserve (χ2 ≥ 12.1, df = 2, P < 0.002). Temporal differences in modeled survival among cohorts were present in both model parameters, b and c (χ2 ≥ 3.9, df = 1, P ≤ 0.04). The minimum and maximum model-predicted survival at age among cohorts within a reserve varied by 1–40%, < 1–27%, and 3–38% at 0.5, 1.0, and 2.5 years of age, respectively. Only in one reserve, Bluff Point, were survivorship curves significantly different among all three cohorts (χ2 ≥ 8.9, df = 2, P ≤ 0.01). Weibull function shape parameters (c) indicated that reserve survivorship schedules followed theoretical type I (e.g., Deep Bay June 2006 cohort), type II (e.g., Ocracoke August 2006 cohort), and type III (e.g., Bluff Point June 2006 cohort) survivorship curves (Appendix 4). No reserve exhibited all three types of survivorship curve.
The Weibull functions were significantly different in space among reserves within a cohort (χ2 ≥ 83.8, df = 2, P < 0.0001; Figure 5). Spatial differences in modeled survivorship were present in b and c (χ2 ≥ 17.6, df = 1, P < 0.0001). Among reserves, cohort survival at age was highly variable, falling in the ranges 48–98%, 36–89%, and 8–53% at 0.5, 1.0, and 2.5 years of age, respectively (Figure 5). Cohorts at Deep Bay typically had the highest probability of survival to a given age (Figure 5). The survival rank order patterns among the remaining reserves varied among cohorts. For instance, Bluff Point, Hatteras, and Crab Hole had the lowest survival probabilities for the June 2006, August 2006, and May 2007 cohorts, respectively (Figure 5).
Density Dependence
In general, oyster growth was density independent in space (i.e., across reserves within a cohort) and time (i.e., across cohorts within a reserve), whereas oyster survival was density dependent in both space and time. When present and significant, density dependence was negative, such that demographic rates decreased with increasing age-0 recruit (survival) or total density (growth).
Growth.—Spatial density dependence in growth was not detectable on trays (Table 4). From quadrat sampling, the strength of spatial density dependence in growth, measured as the slope between mean LVL at age 1 and total density, varied among cohorts (interaction: P = 0.02; Table 4). Only the May 2007 cohort exhibited spatial density dependence at the metapopulation scale (slope: P = 0.005; Figure 6a). The strength of temporal density dependence in growth varied among reserves on both substrates, with total density as the covariate (interaction: P ≤ 0.04; Table 4). Growth was density independent in time at all reserves except Bluff Point and Deep Bay on settlement trays and Ocracoke in quadrats (slope: P ≤ 0.04; Figure 6b, c).
Survival.—Spatial density dependence in survival was not detectable on trays (Table 4). From quadrat sampling, the strength of spatial density dependence in survival, measured as the slope between survival to age 1 and age-0 recruit density, did not vary among cohorts (interaction: P = 0.6; Table 4). With the interaction term removed, survival among pooled reserves significantly decreased as recruit density increased (recruit density: P = 0.004; Table 4), suggesting that both cohorts exhibited spatial density dependence at the metapopulation scale (Figure 7a). The strength of temporal density dependence in survival was similar among reserves (interaction: P ≥ 0.3; Table 4). With the interaction term removed, mean survival differed among reserves on settlement trays (reserve: P = 0.001) but not in quadrats (reserve: P = 0.2; Table 4; Figure 7b, c). More importantly, postrecruitment survival to age 1 was density dependent in time when pooled across cohorts at all reserves for both trays and quadrats (recruit density: P ≤ 0.05; Table 4; Figure 7b, c).
FIGURE 4.
Seasonalized von Bertalanffy growth functions for the (a) June 2006, (b) August 2006, and (c) May 2007 oyster cohorts at the six reserves studied. The percent difference among the reserves in the predicted size at age 2.5 years is shown to the right of each panel. Different superscripts in the legend denote significant differences in the growth functions among the reserves within each cohort. The vertical dashed lines represent the maximum age/date at which the cohorts were observed.
FIGURE 5.
Weibull survivorship functions for the (a) June 2006, (b) August 2006, and (c) May 2007 oyster cohorts at the six reserves studied. See Figure 4 for additional information.
TABLE 4.
Analysis of covariance models testing for temporal and spatial density dependence in oyster growth and survival in tray and quadrat samples. The models included age-0 oyster recruit or total oyster density as covariates and reserve or cohort as categorical factors. P-values are significant at α = 0.05. In settlement trays, recruit density = total density.
DISCUSSION
Over the 3-year study, oyster recruitment and total density at reserves increased 15- and fivefold, respectively. The unprecedented high oyster densities in certain reserves (up to 6,500/m2 at study end) modified demographic rates such that further density increases may be regulated by the prevalent, but relatively weak density-dependent survival. Oyster demographic rates varied significantly among reserves and not in the same manner, such that certain reserves could be classified as the strong “recruiters” (e.g., Crab Hole), others as the fast “growers” (e.g., Ocracoke), and yet others as the high “survivors” (e.g., Deep Bay) based on a relatively consistent (e.g., 2 out of 3 cohorts) top statistical ranking for the respective demographic rate. Superlative performances were altered by cohort dynamics at both intra- and interannual scales, suggesting that the demographic rank order of reserves within the network is a shifting mosaic that may serve as a buffer to biotic and abiotic variability (Gaines et al. 2010).
While identifying the demographically “best” reserves is ideal from a management perspective, our surveys and those conducted by Mroch et al. (in press) suggest that no single reserve concurrently maximized all of the demographic rates over time. Still, the results here can be used by managers in a variety of ways. For instance, the prevalence of densitydependent survival suggests that the oysters in this system are habitat rather than recruitment limited. Consequently, the utility of stock enhancement programs may be minimal and, when implemented, should focus on reserves characterized by high survivorship (e.g., Deep Bay) during times of low recruitment (e.g., fall). The addition of habitat (i.e., artificial reef material) and reserve expansion should focus on reserves such as Bluff Point that are characterized by high recruitment but densitydependent growth and survival. Based on a different suite of demographic rates, Mroch et al. (in press) recommended stock enhancement at Bluff Point and reserve expansion at Ocracoke because these reserves had the highest relative per capita fecundity and reproductive potential per m2 (i.e., integration of per capita fecundity, oyster density, and size structure), respectively. The contradictory management recommendations driven by this apparent demographic mosaic highlight the need to integrate multiple potentially confounding demographic rates within a metapopulation framework to guide comprehensive restoration and management in reserve networks aimed at restoring populations.
FIGURE 6.
Density-dependent growth of oysters in (a) space and (b)—(c) time. Data grouped by reserve within cohorts (panel a) reveal spatial density dependence in growth as a function of total density in quadrats. Data grouped by cohort within reserves reveal temporal density dependence in growth as a function of total density in settlement trays (panel b; age-0 recruit density = total density) and quadrats (panel c). Regression lines and parameters with significant density dependence are shown (see Table 4 for model details).
Density
Our results support the work of others suggesting that the density and age/size structure of protected species can rapidly increase in marine reserves (Halpern and Warner 2002; Halpern 2003; Claudet et al. 2008; Lester et al. 2009; Babcock et al. 2010). We observed a rapid and significant increase in recruit (15-fold), total (fivefold), and legal-sized (threefold) oyster density in 3 years, consistent with the 1–5-year timeframe reported for marine reserves in the literature (Halpern and Warner 2002; Gell and Roberts 2003 and references therein; Claudet et al. 2008; Babcock et al. 2010). Moreover, the recruit and total oyster densities measured at West Bay, Hatteras, and Deep Bay increased 15- to 20-fold over the 5 years since surveys by Powers et al. (2009) and compared favorably with other systems (Figure 8a).
For sessile species such as oysters, recruitment is an important mechanism in determining reserve efficacy and restoration success (Powers et al. 2009). Oyster recruitment, while variable in both space and time, was evident in all years and reserves in our study. The densities of oyster recruits (≤25 mm LVL) in our study (135–1,100/m2; Appendix 2) compared favorably with the recruit densities (≤30 mm LVL) in oyster reserves in Chesapeake Bay (∼350/m2; Schulte et al. 2009). We observed peak oyster recruitment in July—August (May—June cohort), which is consistent with May peaks in per capita fecundity and reserve spawning potential (Mroch et al., in press) as well as the subsequent June peaks in weekly settlement patterns observed in Pamlico Sound (Eggleston et al. 2011). The combination of increased recruitment, rapid growth, and relatively high survival improved reserve age/size structure, which should serve as a positive feedback loop for reserve broodstock function since per capita fecundity increases exponentially with oyster size (Mann and Evans 1998; Mroch et al., in press). Moreover, the prevalence of self-recruitment among reserves (Haase et al. 2012) may complete the positive feedback loop and provide an explanation for the observed rapid increase in oyster recruitment and density over time.
FIGURE 7.
Density-dependent survival in (a) space and (b)—(c) time. Data grouped by reserve within cohorts (panel a) reveal spatial density dependence in survival as a function of age-0 recruit density in quadrats. Data grouped by cohort within reserves reveal temporal density dependence in survival as a function of age-0 recruit density in settlement trays (panel b; age-0 recruit density = total density) and quadrats (panel c). Regression lines and parameters with significant density dependence are shown (see Table 4 for model details).
FIGURE 8.
Panel (a) shows a comparison of mean total oyster density among protected reefs in the York River, Virginia (YR; Nestlerode et al. 2007); the Wicomico River, Virginia (WR; Schulte et al. 2009); Chesapeake Bay tributaries in Virginia (VA; Mann 2010); Back Sound in North Carolina and West Bay, Hatteras, and Deep Bay in Pamlico Sound from 2002 to 2003 (BS and PS*; Powers et al. 2009); and Pamlico Sound (PS; this study). Panels (b) and (c) show a comparison of von Bertalanffy growth functions estimated from this study (shaded area = June 2006 cohort) with historic growth data from Pamlico Sound (Coker 1907) and Chesapeake Bay (Harding et al. 2008) (panel b) and growth functions from oyster reserves in Chesapeake Bay (Paynter et al. 2010), disease resistant and triploid (i.e., sterile) oysters in Chesapeake Bay (Harding 2007), and harvested oysters in Chesapeake Bay (Southworth et al. 2010) (panel c). Panel (d) shows a comparison of Weibull survivorship functions estimated from this study (shaded area = June 2006 cohort) with survivorship data for disease-resistant oysters in Chesapeake Bay (Encomio et al. 2005), triploid oysters in Chesapeake Bay (Calvo et al. 1999), harvested oyster reefs in Chesapeake Bay (Mann et al. 2009; Southworth et al. 2010), and an oyster life history table (Malinowski and Whitlatch 1988).
Growth and Survival
Oyster growth and survival in the reserves surveyed in this study compared favorably with historical and contemporary demographics reported in the literature (Figure 8b, c, d). It should be noted that we measured growth and survival in the upper vertical half of oyster reefs where these demographic rates tend to be highest (Lenihan 1999), so these results may represent demographic potential.
Growth and survival were generally more variable in space (among reserves) than in time (among cohorts within a reserve). For instance, differences in LVL∞ among reserves but not cohorts within a reserve suggest that oyster growth potential depends on where, not when, oyster recruitment occurs. The implications of the observed spatiotemporal variation in demographics are threefold. First, empirical surveys must be conducted at spatiotemporal scales matching demographic variation, which precludes one from applying a universal (in both space and time) age—length key for assessments (sensu Mann et al. 2009; Southworth et al. 2010). Second, the shifting mosaic of reserve demographic rank order highlights a benefit of reserve networks, whereby multiple reserves may be demographically more stable than a single isolated reserve (Gaines et al. 2010). Lastly, the shifting demographic mosaic complicates applied management decisions such as resource allocation and site selection.
Differences in oyster growth and survivorship curves among reserves and among cohorts within a reserve may be caused by several mechanisms (and their interactions), including genetic differentiation, temperature, salinity, competition for resources (e.g., space), disease, and predation (Sebens 1987; Kennedy 1996; Lenihan 1999; Lenihan et al. 1999). Although temperature likely caused seasonal growth stasis (Paynter and DiMichele 1990), tw was invariant across cohorts and reserves, suggesting that temperature is too homogenous soundwide (Roelofs and Bumpus 1953; Durham 2009) relative to oyster tolerance (2–36°C) and optimal conditions (20–30°C; Shumway 1996 and references therein) to have directly caused intra- or interreserve differences in oyster growth and survival. Salinity is likely the overarching driver of demographic variability in oyster growth and survival rates (Wells 1961). Oyster growth rates are typically higher in more saline waters (e.g., Ocracoke and Hatteras; Shumway 1996 and references therein), which explains the stunted oyster growth at Crab Hole, where oysters are exposed to freshets flowing south from Albemarle Sound (Epperly and Ross 1986; Xie and Eggleston 1999; Haase et al. 2012). Salinity indirectly affects survival through the increased diversity and abundance of predators (Wells 1961), which disproportionately prey upon small oysters (White and Wilson 1996). The prevalence of type III survivorship curves at reserves closest to oceanic inlets (e.g., Ocracoke and Hatteras) during peak recruitment and type I survivorship curves at reserves close to low-salinity inputs (e.g., Crab Hole and Deep Bay) supports the positive relationship between predation and salinity.
Density Dependence
We tested for spatial and temporal density dependence in two demographic rates—growth and survival—on standardized (settlement trays) and nonstandardized (quadrat) substrates. Observed oyster densities varied from 6- to 97-fold among reserves within a cohort and from 2- (at Ocracoke and Hatteras) to 34-fold among cohorts within a reserve, providing relatively good contrast for detecting spatial and temporal density dependence (Connell 1985). Density-dependent oyster growth and survival to age 1 was present in space when data were pooled among reserves and time when data were pooled among cohorts, suggesting that density dependence is detectable at annual and metapopulation scales. Temporal density dependence was more prevalent than spatial density dependence, which is significant because temporal, not spatial, density dependence is necessary for population regulation (Stewart-Oaten and Murdoch 1990; Murdoch 1994).
The strength of density-dependent oyster growth (0.006– 0.01) and survival (9 × 10-5 to 1 × 10-4), measured by the magnitude of the slope between recruit or total density and the demographic response, was similar to the values reported by Jenkins et al. (2008) for the acorn barnacle Semibalanus balanoides (growth: 0.001–0.008; survival: 8 × 10-5 to 1 × 10-4) but relatively weak compared with the values reported for reef fishes (survival > 0.003; Schmitt and Holbrook 1999; Steele and Forrester 2005). Within our study, the strength of density dependence was similar between quadrat and settlement trays despite the presence of a preexisting fouling community— barnacles Balanus spp., ribbed mussels Geukensia demissa, macroalgae Codium fragile, and sponges Cliona spp. and Microciona prolifera—in quadrats that may moderate density dependence through mechanisms such as competition for space (Lohse 2002; Boudreaux et al. 2009). Settlement trays, though initially unfouled, provided a standardized substrate that minimized any confounding artifacts of preexisting local conditions and allowed for more precise individual-based censuses, which can be beneficial if the assumptions of quadrat-based modal analysis are not met (e.g., Puckett et al. 2008).
The density threshold below or above which densitydependent growth or survival occurred was not consistent. For example, oysters on settlement trays at Deep Bay exhibited temporal density dependence in growth at total densities 120% lower than those at West Bay, where oyster growth rates were density independent (Figure 6a). Moreover, the strength of temporal density dependence in the survival of oysters in quadrats at Deep Bay and West Bay was equal despite a fourfold difference in recruit density. The inconsistent density-dependent response in reserves harboring disparate oyster densities contrasts with the density-dependent response of oysters in South Carolina (Knights and Walters 2010) and suggests that reserve carrying capacities are highly variable due, most likely, to a suite of site-specific biotic and abiotic mechanisms.
Mechanistically, the significance of recruit and total density in explaining density-dependent survival and growth, respectively, is noteworthy. Recruit density as an explanatory variable suggests that predation is a mechanism underlying density-dependent survival (Schmitt and Holbrook 1999). Mechanical limitations prevent many oyster predators (e.g., mud crabs Panopeus herbstii) from consuming large oysters; thus, predatory responses are largely determined by recruit (i.e., prey) density rather than total oyster density (Rindone and Eggleston 2011). Total oyster density as an explanatory variable suggests that competition for space is a mechanism underlying densitydependent growth (Sánchez Lizaso et al. 2000). Because oysters are sessile organisms, both intra- and intercohort density increases physical contact which, in plate tectonic fashion, can lead to subduction of some oysters and subsequent reductions in growth rate (authors' personal observations).
Conclusions
This study quantified spatiotemporal variation in key oyster demographic rates in a no-take reserve network that, when coupled with evidence of larval connectivity among reserves in the network (Haase et al. 2012), provides support for applying a metapopulation conceptual framework in this study system. From an oyster restoration perspective, our work supports the recent work by Schulte et al. (2009), Powers et al. (2009), and Paynter et al. (2010) indicating that establishing no-take reserves is an effective and ecosystem-based approach for developing a protected broodstock as part of an oyster restoration portfolio. Yet, a one-size-fits-all restoration strategy is likely to be misguided. The spatiotemporal variation in demographic rates and density dependence are both important when one is considering a potential restoration strategy due to their importance in determining species biomass within reserves and the broodstock function of reserves.
ACKNOWLEDGMENTS
We are especially grateful to the following for their help with field collections: R. Mroch, R. Rindone, G. Plata, A. Haase, G. Bell, C. Durham, E. Millstein, J. Wiggs, M. Moorman, M. Radlinski, and J. Eggleston. We are also grateful to the North Carolina Division of Marine Fisheries for their logistical support (C. Hardy and S. Slade), as well as to J. Taylor for image analysis software training, G. Bell for statistical guidance, and K. Shertzer, N. Haddad, and two anonymous reviewers for critically reviewing an earlier version of this manuscript. Funding for this project was provided by a NMFS Sea Grant population dynamics graduate fellowship (to B. Puckett; NA08OAR4170764), a North Carolina Sea Grant (to D. Eggleston), the North Carolina Fisheries Resource Grant Program (to D. Eggleston; 06-EP-03, 07-EP-04), the Blue Crab Advanced Research Consortium via the University of Maryland (to D. Eggleston), an NSA Michael Castagna Student Grant for Applied Research (to B. Puckett), and a Raleigh Saltwater Sportfishing Club Dale Ward Memorial Scholarship (to B. Puckett).
REFERENCES
Appendices
APPENDIX 1: LENGTH FREQUENCY HISTOGRAMS
FIGURE A. 1.1. Length frequency histograms obtained from n replicate quadrat samples at West Bay, Ocracoke, Hatteras, Crab Hole, Bluff Point, and Deep Bay oyster reserves from June 2006 to July 2008. The best-fitting modes (dashed lines) and cumulative modes (solid lines) are overlaid. The June 2006 (J06), August 2006 (A06), and May 2007 (M07) cohorts are identified in the panels. See Table 2 in the text for cohort modal means and fit statistics. Note the differences in y-axis scales.
APPENDIX 2: OYSTER DENSITIES
TABLE A.2.1.
Mean ± SE total, recruit (≤25 mm LVL), and legal (>75 mm LVL) density at six oyster reserves from June 2006 to July 2008. Different superscripted letters and numbers denote significant differences in mean density among reserves within a given sampling date and among sampling dates within a given reserve (i.e., simple effects), respectively.
APPENDIX 3: GROWTH FUNCTIONS
FIGURE A.3.1. Observed left valve length (LVL) at age and date (open circles) of three cohorts—June 2006, August 2006, and May 2007—of oysters obtained from settlement trays at six reserves. The solid lines represent the best-fitting seasonalized von Bertalanffy growth functions. Parameters are as follows: LVL∞ is the maximum asymptotic LVL for the ages sampled, K is a curvature parameter that determines how fast LVL∞ is approached, t0 is the theoretical age (year) when LVL is zero, C is related to the magnitude of the seasonal oscillation, and tw is the time at which growth is lowest. Asterisks indicate that parameters were constrained to reported values. The dashed lines represent the upper and lower 95% confidence limits. The dotted lines represent the ages at which cohorts are predicted to attain the commercial size limit (76 mm). Modal means (closed squares; ± SEs) obtained from length frequency analysis of quadrat samples are overlaid.
APPENDIX 4: PROPORTION SURVIVING FUNCTIONS
FIGURE A.4.1. Proportion surviving to age and date (open circles) of three cohorts—June 2006, August 2006, and May 2007—of oysters obtained from settlement trays at six reserves. The solid lines represent the best-fitting Weibull survivorship functions. Parameters are as follows: b is the scale parameter and c is the shape parameter; values of c > 1, 1, and <1 correspond to type I, II, and III survivorship curves, respectively. The dashed lines represent the upper and lower 95% confidence limits. The dotted lines represent the ages at which 50% of the cohorts are predicted to survive. The proportions surviving to ∼age 1 (closed squares; ± SEs) obtained from length frequency analysis of quadrat samples are overlaid for the June 2006 and May 2007 cohorts. Quadrat sampling was not conducted in October, so estimating the August 2006 initial cohort density (and thus survival) was not possible.
