Effect of oyster aquaculture on seagrass Zostera marina at the estuarine landscape scale in Willapa Bay , Washington ( USA )

Both seagrasses and bivalve shellfish provide valuable ecosystem services in estuaries worldwide. Seagrasses are protected by no-net-loss provisions in US federal and state regulations, resulting in precautionary management that avoids any direct impacts from development activity, including shellfish aquaculture. Recent research suggests that oyster aquaculture has direct impacts on native seagrass (eelgrass Zostera marina) at small spatial and short temporal scales in US west coast estuaries. We quantified impacts of oyster aquaculture on Z. marina at the estuarine landscape scale in Willapa Bay, Washington. A model of Z. marina cover outside of aquaculture was created using distance to estuary mouth, distance to nearest channel, salinity, elevation, and cumulative wave stress as factors, and was then used to predict Z. marina distribution within oyster aquaculture beds and compared to an inverse distance interpolation of points outside of aquaculture. The amount of Z. marina cover observed within oyster aquaculture beds was less than predicted, but represented <1.5% of the total predicted amount of Z. marina cover in Willapa Bay in any year. Type of oyster culture bed did not contribute to observed variation, but mechanically harvested beds had significantly less Z. marina cover than beds harvested by other methods. The majority of beds had 65−145% of the model-predicted Z. marina cover and exhibited relatively low variability between years, suggesting that Z. marina as habitat is resilient to oyster aquaculture as a disturbance and does not result in persistent effects at the landscape scale in this estuary.


INTRODUCTION
Both seagrasses and reef-forming bivalves serve a variety of important ecological functions in estuaries, including: enhancing biodiversity and providing structured nursery habitat and refuge from predators for fish and invertebrates (e.g.seagrasses: Bostrom et al. 2006, Gillanders 2007, Heck et al. 2008;bivalves: Coen et al. 1999, Gutierrez et al. 2003, Grabowski et al. 2008), water column filtration and water property enhancement (particularly bivalves as phytoplankton grazers; Prins et al. 1998, Newell 2004, Ferreira et al. 2011), sediment accretion and erosion control through current modification (seagrasses: Peralta et al. 2008, Koch et al. 2009;oysters: Smith et al. 2009, Scyphers et al. 2011), carbon sequestration (seagrasses: Mateo et al. 2007, Fourqurean et al. 2012), and finally, as foraging areas for waterfowl and shorebirds (both seagrasses and bivalves; Caldow et al. 2007, Rivers & Short 2007, Anderson & Lovvorn 2012).Reef-forming bivalves and seagrasses are also declining worldwide due to numerous anthropogenic disturbances, including eutrophication, wetland filling and diking, fishing, and dredging (Orth et al. 2006, Waycott et al. 2009, Beck et al. 2011, Zu Ermgassen et al. 2012).These declines have led to efforts to protect and restore seagrass and bivalve populations, as well as an increased focus on identi-ABSTRACT: Both seagrasses and bivalve shellfish provide valuable ecosystem services in estuaries worldwide.Seagrasses are protected by no-net-loss provisions in US federal and state regulations, resulting in precautionary management that avoids any direct impacts from development activity, including shellfish aquaculture.Recent research suggests that oyster aquaculture has direct impacts on native seagrass (eelgrass Zostera marina) at small spatial and short temporal scales in US west coast estuaries.We quantified impacts of oyster aquaculture on Z. marina at the estuarine landscape scale in Willapa Bay, Washington.A model of Z. marina cover outside of aquaculture was created using distance to estuary mouth, distance to nearest channel, salinity, elevation, and cumulative wave stress as factors, and was then used to predict Z. marina distribution within oyster aquaculture beds and compared to an inverse distance interpolation of points outside of aquaculture.The amount of Z. marina cover observed within oyster aquaculture beds was less than predicted, but represented <1.5% of the total predicted amount of Z. marina cover in Willapa Bay in any year.Type of oyster culture bed did not contribute to observed variation, but mechanically harvested beds had significantly less Z. marina cover than beds harvested by other methods.The majority of beds had 65−145% of the model-predicted Z. marina cover and exhibited relatively low variability between years, suggesting that Z. marina as habitat is resilient to oyster aquaculture as a disturbance and does not result in persistent effects at the landscape scale in this estuary.KEY WORDS: GIS • Eelgrass • Aerial photography • Shellfish aquaculture • Estuary OPEN PEN ACCESS CCESS fying and mitigating the factors that are negatively impacting each of them (Brumbaugh & Coen 2009, Schulte et al. 2009, Orth et al. 2010a, Thom et al. 2012).
Bivalve shellfish aquaculture acts as an anthropogenic disturbance to seagrass, but can also positively interact with seagrass and may restore some of the services lost where native bivalve populations have declined (reviewed by Dumbauld et al. 2009, Forrest et al. 2009, Coen et al. 2011, McKindsey et al. 2011).The interactions between bivalve shellfish aquaculture and seagrasses have been widely studied at the experimental scale, particularly those between Pacific oysters Crassostrea gigas, cultured in many estuaries worldwide, and eelgrass Zostera marina, a common temperate species of seagrass.Shellfish influence seagrass via 3 primary mechanisms: physical structure of shells/reefs, nutrient addition to sediments and water column via excretion, and increased water clarity via filtration.Oysters are expected to compete for space with seagrass, especially when they are cultured directly on the sediment surface.Research on oyster aquaculture in Willapa Bay, Washington on the US west coast suggests that this interaction is nonlinear and thresholds occur above which shoot density of Z. marina declines markedly (Wagner et al. 2012).Oysters grown on structures can have additional physical effects, including shading and sediment erosion around the structure.These factors caused 75% reduction in eelgrass cover relative to controls for oyster stake culture and up to 100% loss of eelgrass under oyster racks (Everett et al. 1995), yet oysters grown on longlines with more open space caused little reduction in eelgrass density and cover (Wisehart et al. 2007, Tallis et al. 2009).Further study suggested that eelgrass metrics scaled with spacing between these oyster lines and that both shading and dessication from stranding over the lines contributed to the effect (Rumrill & Poulton 2004).Experimental evaluations of the effects of shading due to oyster culture in suspended bags and hanging culture of floating oyster bags showed reduced eelgrass structure, morphometrics and photosynthesis observed at 26% reduction in subsurface irradiance, but design of the structures was clearly important, as the floating bags only reduced eelgrass directly below the structures (Bulmer et al. 2012, Skinner et al. 2014).Bivalves have also been shown to enhance seagrass growth by supplying nutrients via biodeposits and by improving water clarity via filtration, especially where nutrients in sediment limit seagrass growth and in eutrophic waters where phytoplankton blooms cause shading effects (Peterson & Heck 2001, Booth & Heck 2009, Burkholder & Shumway 2011).Though nutrients in sediment porewater were enhanced by the presence of oysters in Willapa Bay (Wagner et al. 2012) and oysters measurably cleared the water (Wheat & Ruesink 2013), only eelgrass shoot size was affected by oysters at the local scale (Wagner et al. 2012), suggesting that seagrass response differs depending on estuarine conditions.
Seagrasses have been shown to be sensitive to a wide variety of pulse disturbances with parallels to mechanical implements used to harvest shellfish (e.g.boat propellors, anchors, and moorage chains: Dawes et al. 1997, Thom et al. 1998; dredge and fill operations and simple trampling: Erftemeijer & Lewis 2006).Shellfish harvest practices have been less studied, but mechanical harvest implements directly removed plants and generally caused more disturbance than hand harvest or off-bottom longline oyster culture techniques in Willapa Bay (Wisehart et al. 2007, Tallis et al. 2009).
Most of the experimental studies outlined above describe effects of oyster aquaculture on seagrass at small spatial and short temporal scales.There has been an extensive amount of effort devoted to using remote sensing to understand seagrass dynamics at larger estuarine landscape scales, particularly in estuaries where other anthropogenic disturbances such as eutrophication are more likely to affect change (Kendrick et al. 2000, Dekker et al. 2007, Orth et al. 2010b, Lyons et al. 2013), but only recently have investigators addressed impacts of shellfish aquaculture on seagrass at this scale (Ward et al. 2003, Carswell et al. 2006, Barille et al. 2010, Martin et al. 2010, Bulmer et al. 2012).Willapa Bay provides a unique opportunity to examine this interaction at the estuarine landscape scale because shellfish aquaculture within eelgrass was not restricted prior to 2007 when a new permit (US Nationwide Permit 48) was issued by the US Army Corps of Engineers (US ACOE) with protection for native eelgrass Z. marina.This permit did not influence the placement or activities on privately owned or leased aquaculture beds, but if these activities occurred in eelgrass, 'pre-construction notification' was required and any expansion of aqua culture must leave buffers for eelgrass.Initial landscape-scale estimates from Willapa Bay clearly de picted reductions of eelgrass cover on individual beds, but suggested that when averaged over space at similar tidal elevations, the proportion of area with eelgrass present inside and outside of oyster culture grounds was similar (Dumbauld et al. 2009).Here we further that effort for this estuary by first constructing a model to predict seagrass distribution outside of aquaculture using several factors that we suspected could influence the distribution of Z. marina and for which we had spatial data to create layers, including distance to estuary mouth, distance to nearest channel, salinity, elevation, and cumulative wave stress.We then use these factors to predict Z. marina distribution for each aquaculture bed and compare the model-predicted, interpolated, and actual quantities of Z. marina.Our study had several related objectives: (1) quantify the distribution of Z. marina and oyster aquaculture in this estuary, (2) quantify the overall impact of oyster aquaculture on Z. marina, (3) determine the relative impact of different oyster aqua culture harvest methods and bed types, and (4) determine whether any impacts of oyster aquaculture on Z. marina were chronic or transitory by analyzing data from 3 separate years.

Study site
Willapa Bay, Washington (46°N, 124°W; Fig. 1) is a macrotidal estuary on the US Pacific coast, with a large tidal exchange and a well-mixed water column heavily influenced by tidal exchange near the estuary mouth, especially during summer months (Banas et al. 2007).It is the third largest estuary on the US west coast, encompassing 358 km 2 , with 227 km 2 or 62% of the bay being intertidal (Hedgepeth & Obrebski 1981, see Table 1).Nine small rivers contribute to Willapa Bay, with a total watershed area of 2776 km 2 .Shellfish aquaculture occurs on both leased tidelands and privately owned tidelands in US west coast estuaries, with up to 22% of the intertidal estuarine tidelands being devoted to this activity in Willapa Bay (Feldman et al. 2000, Dumbauld et al. 2009, 2011).Non-native Pacific oysters Crassostrea gigas are the primary species cultured, and this single estuary produces approximately 17% (2585 t) of the annual US oyster product (NMFS 2014, M. Morningstar, Washington Dept.Fish and Wildlife, pers.comm., K. Ra mey, California Dept.Fish and Game, pers.comm.).Although seagrasses are in decline worldwide, several US Pacific coast estuaries are charac terized by expansive meadows of seagrass.
There are a variety of techniques used for farming C. gigas in Willapa Bay, including longlines, stakes, and bags, but the vast majority of oysters are grown directly on-bottom (Dumbauld et al. 2011).Three main techniques are used for setting and growing very young oysters: setting larvae on bags of shell in a hatchery, capturing wild spawned larvae on bags of shell, and capturing wild spawned larvae on loose shell distributed on the tideflat.These new oysters are placed on 'seed' beds, usually located further from the mouth of the bay and designed to minimize mortality and loss.After obtaining a sufficient size, they are often transferred to 'fattening' beds, usually located closer to the mouth and designed to maximize growth and meat production.After reaching marketable size, oysters are harvested either by hand-picking them into bushel baskets and large tubs at low tide or by mechanical harvest (metal dredges) when the beds are submerged.Pacific oysters harvested for the fresh shucked market are generally ready for harvest 3−4 yr after planting.

Eelgrass
We conducted an on-ground survey of 4238 stations in Willapa Bay in 2006−2007, where we collected habitat information, including Z. marina cover.Sampling was conducted on a grid, with stations located at 200 m intervals across the intertidal, and each accessible location was visited at low tide by hovercraft.At each station, we recorded the presence and density of 2 seagrass species (native Z. marina and introduced Z. japonica), macroalgae, live oysters and shell in a 10 m 2 area.Density of each eelgrass species, macroalgae (classified loosely into brown, green, red), oyster and shell was classified into 4 categories: absent, present−minor (<25% cover), present−medium (25−75% cover), and present− major (> 75% cover).A pilot survey suggested these categories related directly to quantified assessments using smaller quadrats, but the larger 100 m 2 scale was most useful to classify vegetation signatures in aerial photographs.Additional data on burrows within a 0.25 m 2 quadrat and burrow occupants (clams, polychaetes, and thalassinid burrowing shrimp), as well as basic sediment characteristics, were also taken.Data were recorded into a Trimble GeoXT ® mapping-grade GPS system.At each station, we recorded elevation relative to height above ellipsoid using an L1 Trimble ProXR ® mounted at a fixed position on the hovercraft to create a complete digital elevation model by combining our low-resolution elevation model with a LiDAR dataset previously collected for the upper intertidal.These survey data were post processed, and exported as an ESRI shapefile using Trimble GPS Pathfinder Office v.4.20.We determined the cover and distribution of Z. marina for all areas of the bay by: (1) collecting, georeferencing, and orthorectifying 0.25 m color infrared aerial photography for all of Willapa Bay, WA (1:20 000 with 40% overlap) in May 2005 (USDA-ARS), July 2006 (WA DNR), and June 2009 (USDA-ARS); (2) extracting the mean values from each color band in a 5 m radius buffer at each of the 4238 stations surveyed in 2006− 2007; (3) exploring and selecting the best model to predict the relationship between on-ground estimates of presence (absent, present−minor, present− medium, and present−major) and the best performing relationship between bands (normalized difference vegetation index [NDVI], infrared-red, and infrared-green); and (4) using this model (infrared-red performed best) to predict the probable cover of Z. marina for each pixel in the imagery (mean of 1 × 1 m values).All data were imported into R (R Development Core Team 2013) for analysis, maps were generated using the rgdal (Keitt et al. 2012) and raster (Hijmans & van Etten 2012) packages, and geostatistics calculated using the rgeos and gstat packages (Pebesma 2004, Bivand & Rundel 2014).The photo-extracted cover layers were used as the actual cover and distribution of Z. marina for 2005, 2006, and 2009.Descriptive distributional statistics were then calculated by intersecting the raster of Z. marina cover (0 to 1), with binned ranges of the 5 factor rasters (see 'Other factors' below for description).Summing the raster val-ues resulting from this intersection gives the total amount of Z. marina cover predicted in each band (i.e.Z. marina between 0.3−0.6 m mean lower low water [MLLW]; see Fig. S1 in the Supplement at www.int-res.com/articles/ suppl/ q007p029_ supp.pdf for a more detailed flowchart).
As the elevation increases above ~+1.2m MLLW, several vegetative patterns occur that make predicting Z. marina less accurate: (1) Z. japonica, an introduced species of seagrass, begins to appear.This seagrass can be very dense, but is difficult to quantify using aerial infrared photography due to its short and narrow profile.However, at high densities, it has a signature similar to low-density Z. marina.(2) Z. marina can occur in pools at upper elevations and in narrow drainage channels that are not detectable in the elevation layer.(3) In 2005 and 2006, there was an active removal program in place for Spartina alterniflora, an introduced cordgrass that occurs in the upper intertidal (only 5% occurs below +1.2 m MLLW, our survey results; Civille et al. 2005).Herbicide-treated S. alterniflora has a similar signature as medium-density Z. marina, and unsprayed S. alter niflora has a similar signature to high-density Z. marina.Due to these challenges of accurately identifying Z. marina at upper elevations, and be cause the vast majority of active oyster aquaculture occurs below +1.2 m MLLW, we chose to focus on analyzing the interactions of aquaculture and Z. marina below this tidal elevation.

Aquaculture
We conducted interviews with the major clam and oyster shellfish growers in Willapa Bay in 2006 and 2007.They provided location, species, aquaculture bed type, and harvest method information for each aquaculture bed (n = 458).This information was entered into a GIS.For this analysis, we limited exploration to the 282 active on-ground oyster beds with > 90% of the bed in the intertidal (above -1.5 m MLLW) and 95% of the bed below +1.2 m MLLW.Descriptive statistics were calculated for oyster aquaculture distribution using each of the 5 predictor layers.
A 5 × 5 m raster of oyster aquaculture (0/1 representing absence/presence) was created.Multiple temporary rasters were created for each of the 5 factor rasters (see 'Other factors' below for description) representing binned ranges of values for each factor.For example, 17 temporary rasters were created for each 0.3 m elevation interval between −1.5 and + 3.7 m MLLW, with cells coded as 1 if they were in this interval and 0 otherwise.The aquaculture raster was multiplied by each temporary raster.Summing the raster resulting from this intersection (counting the number of cells with a value of 1) gives the total area of aquaculture in each band (see Fig. S2 in the Supplement for a more detailed flowchart).

Other factors
We chose 5 abiotic factors that we suspected could influence the distribution of Z. marina and for which we had spatial data to create layers.We created 5 × 5 m resolution rasters for each of these factors, including: (1) elevation, reported in meters and feet relative to MLLW; (2) distance to the estuary mouth, in meters using cumulative cost distance; (3) distance to the nearest channel, in meters using Euclidean distance to the nearest channel with water present in aerial photos; (4) cumulative wave stress, a sum of the 6 min wind speeds over a year using direction, fetch distance, elevation, and tide stage to calculate a relative exposure to wind-driven wave scour; and ( 5) salinity (ppt) using a bay-wide prediction from a model of the 5th quantile wet season salinity levels from 5 long-term monitoring locations by distance to estuary mouth (www.ecy.wa.gov/programs/eap/ mar_ wat/ data.html).These factors had variable resolution but were standardized to allow prediction across the entire bay.The elevation raster was created by interpolating high-resolution LiDAR elevation data (http:// pugetsoundlidar.ess.washington.edu/) for the upper intertidal with 200 m grid elevation data from the lower intertidal (USDA-ARS) into a 10 × 10 m raster (www.onrc.washington.edu/GIS/MDAIportal.html).This raster was resampled to a 5 × 5 m grid for this analysis.

Model and predicted Z. marina layers
We used 2 methods to create species distribution models (SDMs) and predict Z. marina levels in oyster aquaculture areas.We distributed ~5000 random points across the intertidal of Willapa Bay between −1.5 and +1.2 m MLLW and outside of any type of aquaculture, including active and inactive beds, clam beds, longlines, and seed racks.At each of these locations, the actual cover of Z. marina in 2005, 2006, and 2009 was extracted and averaged from a 5 m radius buffer around each point.
For the first method, the predicted Z. marina layer was created by interpolating the actual Z. marina probabilities at each of the random points.Each pixel in this layer represented the inverse-distance weighted mean of the 10 nearest random points.Therefore, pixels near points with high Z.marina densities would be predicted to have high Z.marina densities as well.
For the second method, the values from the 5 predictor layers (elevation, distance to mouth, distance to nearest channel, wave stress, and salinity) were also extracted from the random points and used to generate a logistic generalized additive model (GAM) of the relationship between the predictor layers and Z. marina cover for each year (mgcv package in R; Wood 2006Wood , 2011)).The best model for each year was selected using the root mean square error (RMSE) of a 10-fold cross-validation.We expected that this model would provide a more reliable and spatially precise measure of eelgrass presence if the factors represented by predictor layers were responsible for eelgrass distribution.The selected model and the predictor layers were then used to predict Z. marina across the entire tideflat -including active on-ground oyster aquaculture areas -for each year.The value of each cell represents the probability of observing Z. marina in this pixel area, i.e. a pixel value of 0.10 represents a 10% probability of Z. marina occurrence in that pixel.

Extractions and calculations
Active on-ground oyster aquaculture beds were overlain on the actual, interpolated, and modelpredicted Z. marina layers for each year, and the sum of the pixel values were extracted for each aquaculture bed -giving the total quantity of Z. marina actually observed, interpolated, and model-predicted for each bed for each year.These values were attached to the corresponding bed in the oyster aquaculture table for analysis.

Estuary totals
To calculate the total effect of oyster aquaculture on Z. marina, we took the sum of the actual, predicted, and interpolated probabilities for all of the oyster beds for 2005, 2006, and 2009.Subtracting the actual amount of cover from the predicted and interpolated amounts gives an estimate of the total amount of eelgrass cover missing due to oyster aquaculture.To visualize results across oyster bed size, we plotted log-transformed amounts of observed Z. marina cover versus log-transformed predicted and interpolated cover amounts and fit a linear relationship with 95% confidence interval and 95% prediction interval for each year.Oyster beds below the 1:1 line have less Z. marina cover than predicted, while beds above the 1:1 line have more Z.marina cover than predicted.

Oyster culture bed type and harvest method
To control for differences in oyster bed size, we calculated the proportion of Z. marina predicted or interpolated that was actually observed for each bed.We fit mixed effects models (nmle package in R; Pinheiro et al. 2014) with main effects of bed type and harvest method and random effects of bed and year (to control for pseudo-replication) to this proportion.Pairwise comparisons (multcomp package in R; Hothorn et al. 2014) were made for significant main effects and models were compared using the Bayesian information criterion (BIC).

Temporal patterns
Using the proportional values calculated above, we calculated the mean and standard deviation of the proportion of actual observed to model-predicted Z. marina for each oyster bed across each of the 3 sample years.We predicted that mechanically harvested beds would either exhibit chronically low proportions of Z. marina, if the effects of dredging are long-lived, or high variability, due to a rapid removal (mechanical harvest) and recovery (regrowth), relative to more stable hand-picked beds.

General distribution of seagrass and shellfish aquaculture
About 55% of the intertidal area in Willapa Bay is under private ownership, and our data indicates that bivalve aquaculture occupied 23% of the intertidal area (Table 1).This included both clam and oyster aquaculture beds, seed racks, and oyster seed beds that were only partially or occasionally used.Active intertidal oyster aquaculture occupied just 8% of the intertidal.We observed a vegetation signal in aerial photographs on 41% of the intertidal (between −1.5 and + 3 m MLLW).Of this, 80% occurred below +1.2 m MLLW and was predominately Zostera marina.
The distribution of oyster aquaculture at lower tidal elevations clearly overlapped that of Z. marina (Fig. 2a).While Z. marina was broadly distributed over a large range of distances from the estuary mouth (Fig. 2b), active oyster aquaculture was largely restricted to the central portion of the bay.A large portion of the south end of the bay (farther from the mouth) is partially or opportunistically used for growing juvenile oyster seed (seed beds, but not taken into account in Fig. 2).Distributional patterns of tideflat area, active oyster aquaculture, and Z. marina by other predictor variables (salinity, distance to nearest channel, and stress) were skewed, with most aquaculture and eelgrass being present at higher salinity (>12 ppt), close to the channel, and at relatively low cumulative wave stress (Figs.S3−S5 in the Supplement at www.int-res.com/articles/ suppl/ q007 p029_ supp.pdf).

Model and predicted aquaculture effects
All of the predictor variables contributed significantly to describing the distribution of Z. marina in a logistic GAM (Table 2).Models including elevation, salinity, and distance to estuary mouth as covariates and stress and distance to nearest channel as independent effects were found to be the best fit for years 2005, 2006, and 2009 (10- The number of aquaculture beds with more Z.marina present than predicted equaled that with less eelgrass present than the model predicted in all 3 years (2005 shown in Fig. 3a; 2006 and 2009 in Fig. S6 in the Supplement; both axes log-transformed).A group of oyster beds, however, fell well below and outside the 95% confidence and prediction intervals around a 1:1 line, which represented beds that had the same amount of Z. marina cover as predicted, perhaps representing chronic effects (see mixed effects model analyses below).Results were similar for interpolated values (Fig. 3b).The total area of Z. marina estimated to be missing using a model prediction in 2005 and 2006 was only 22 and 8 ha, respectively (Table 3).In 2009, there were 0.4 ha, more Z.marina present than predicted by the model.The total area of Z. marina estimated to be missing using the interpolation prediction was higher for all years, at 80, 84, and 60 ha, respectively (Table 3).Although large in aggregate, even the highest estimate is <1.5% of the total amount of Z. marina cover found in Willapa Bay in these 3 years.

Aquaculture effects by harvest method and oyster culture bed type
Oyster bed type (direct, fattening, longline, seed, and mixed) was not a significant contributor to the variation in the proportion of Z. marina actually observed versus model-predicted and interpolation-predicted eelgrass present (lower BIC value for a mixed effects model without bed type; Table 4).Harvest method (hand, mechanical, mixed), however ,was found to be significant (ANOVA, p < 0.001), and a Tukey's post hoc test revealed that mechanically harvested beds have a significantly lower amount of Z. marina cover relative to beds that are harvested by hand or where a mixed harvest technique is employed (Table 5).
Mechanical harvested beds on average had 100% and 92% of the modeland interpolation-predicted Z. marina, respectively, hand-harvested beds had 120% and 127% of the predicted and interpolation-predicted Z. marina, and mixed harvest beds had 117% and 97%, respectively.Trends for some individual oyster culture beds were quite evident in aerial photographs (Fig. 4), and when examined over time, exhibited 4 general patterns: (A) beds with chronically low Z. marina and low variation across years, (B) beds with the expected levels of Z. marina with low variation across years, (C) beds with highly variable levels of Z. marina across years, and (D) beds with high levels of Z. marina with high variation across years (Fig. 5).The majority of beds exhibited expected levels of Z. marina with low variation across years.All of the beds with < 65% of the mean expected amount of Z. marina cover present (n = 24) were mechanically harvested beds and demonstrated a chronically low level of Z. marina cover present across years.Most of the beds displaying high variability were handharvested or those where mixed harvest techniques were used.

Overall effects
Our results demonstrate a negative effect of oyster aquaculture on the native seagrass Zostera marina at the landscape scale in Willapa Bay, WA, but also show that this impact is small compared to the overall signature of both Z. marina and oyster aquaculture in this estuary.For its relatively small size, Willapa Bay has a proportionally large intertidal area (62% exposed on extreme low tides, above −1.5 m MLLW) compared with other estuaries where shellfish aquaculture is important on the US west coast and elsewhere in the world (Dumbauld et al. 2009).Oyster aquaculture occurs on a relatively large portion of this tideflat (8%, or 1764 ha; Table 1), and we found Z. marina was also important, covering 38% of the tideflat (6951 ha).These estuarine areas are also proportionally large compared to other estuaries, and the tidal elevation where oyster aquaculture and Z. marina are found overlap substantially (Fig. 2a), suggesting the potential for significant interaction.The amount of eelgrass cover observed within oyster aquaculture beds was generally less than the amount predicted by a logistic GAM built to describe the distribution of presumably undisturbed eelgrass outside of aqua-  culture.This estimated impact of aquaculture on Z. marina was slightly greater when compared to a simple interpolation of values from random points outside of aquaculture.However, while the total area of Z. marina declined slightly over time in our study, <1.5% of either the total predicted or interpolated amount of Z. marina cover was missing (maximum of 80 ha) and could thus potentially be attributed to aquaculture in any single year.This lack of substantial overall impact is similar to the few studies conducted at the estuarine landscape scale elsewhere, but these studies focused only on effects outside of the direct aquaculture signature.Ward et al. (2003) for example found no detectable effect of rack-cultured oysters (137 ha) on Z. marina (2390 ha) in Bahia San Quentin, Mexico using maps created from satellite imagery.Similarly, no observable effect of hanging-bag oyster culture was observed using remote sensing in New Zealand where Z. muelleri expanded outside of the direct signature of the lines (5 ha of lines, eelgrass increased from 610 m 2 to 2537 m 2 ; Bul-mer et al. 2012).Results are also consistent with mapping efforts in France, where traditional onbottom and hanging-bag oyster culture had no observable effects on Z. noltii or Z. marina outside the direct signature of aquaculture (i.e.below the structures; Barille et al. 2010, Martin et al. 2010).Our data suggests that the effect is relatively small even within the direct aquaculture signature.Though use of several models at the same time to distinguish effects has been recommended, GAMs have been successfully used to create SDMs and predict seagrass distribution, and they tend to give a large precautionary estimate of presence relative to other methods (Kelly et al. 2001, Bekkby et al. 2008, Downie et al. 2013, Schubert et al. 2015).Our models only explained about 50% of the variation in eelgrass distribution.Models with ele vation, salinity, and distance to estuary mouth as covariates and stress and distance to nearest channel as independent effects were found to produce the best fit.Some of these predictors, such as tidal elevation and distance to estuary mouth, which were the most significant, were likely only proxies for mechanistic variables such as light and dessication, which have been shown to limit the lower and upper distribution of these plants in estuaries re spectively (Dennison 1987, Koch 2001, Boese et al. 2005, Thom et al. 2008), but for which we did not have spatial data.While the predictor variables used in the model are also auto-correlated, our goal here was to compare the distri bution of Z. marina within and outside oyster aquaculture, so the assumption of independence is less important and perhaps interaction terms be tween these factors and aquaculture are of primary interest.

Oyster culture bed type
We found that the type of oyster culture bed (direct, fattening, longline, seed, and mixed) did not explain any of the variation observed in the proportion of Table 5. Coefficients from mixed effect models using restricted log-likelihood and Tukey's post hoc test (small letters in parentheses) to compare harvest methods.Dredge harvest beds had significantly lower proportions of actual to model and interpolation-predicted Zostera marina than hand-harvested beds.However, dredge harvest beds had 100% of the model-predicted and 92% of the interpolation-predicted Z. marina, while hand-harvested beds had 120% and 127% of the model-and interpolation-predicted amount of Z. marina cover, respectively Z. marina observed relative to the interpolationpredicted or model-predicted amount of cover at the landscape scale.Oysters are typically raised in a hatchery and planted as seed (small spat on cultch shells) either directly on the bottom or suspended on longlines in Willapa Bay.Some growers plant seed on beds that are strictly used for seed during the first year of growth and then transfer these oysters to fattening beds with better growing conditions where they grow for a final period before harvest.Seed beds are usually located in the upper portion of the estuary, while fattening beds are located closer to the estuary mouth.Some oyster culture beds known as direct harvest beds are planted with seed and oysters allowed to grow for the full 3−4 yr cycle or until harvest.The amount of oyster cover thus increases on these beds as they grow and also depends on initial planting density, which varies by location and grower preference (typically 10−20% cover, or 250−750 cultch bags ha −1 ).Wagner et al. (2012) showed that competition for space between oysters and seagrass did not follow a 1:1 relationship and a threshold of about 20% shell cover existed, above which Z. marina density declined exponentially, at least at small experimental scales.With the exception of areas where oysters are allowed to form reefs or hummocks near the south end of Willapa Bay, however, cultured oyster density rarely reaches this threshold (Dumbauld et al. 2006, 2011, Tallis et al. 2009), and eelgrass commonly co-exists with oysters throughout the growout cycle, albeit at reduced density.Despite measurably enhanced sediment or ganic concentration and porewater nutrients due to the presence of oysters in Willapa Bay (Tallis et al. 2009, Wagner et al. 2012), and the potential positive effects this has been shown to have on seagrasses elsewhere (Peterson & Heck 2001, Carroll et al. 2008, Booth & Heck 2009), Z. marina shoot growth was shown to be little affected by the presence of oysters in Willapa Bay experimental studies (Tallis et al. 2009, Wagner et al. 2012).We suspect that on balance, the effect of bottom-cultured oysters on eelgrass in Willapa Bay was variable enough at smaller spatial scales to eliminate any significant effect at the larger landscape scale in our study.Oysters grown on structures would not be expected to have the direct physical effects that placement on the bottom has on seagrass, but instead have effects such as shading and sediment erosion beneath and around the structures themselves.Oysters grown on longlines represented only 4% of the culture area in our study, and again we observed no difference relative to the model-predicted or interpolation-predicted amount of eelgrass cover in these areas at the landscape scale.Effects clearly depend on the structure, with densely planted stakes and oysters grown in bags on racks substantially reducing eelgrass cover in a previous US west coast study (75% and 100% loss, respectively, relative to controls; Everett et al. 1995), yet oysters grown on longlines with more open space caused little reduction in eelgrass density and cover (Wisehart et al. 2007, Tallis et al. 2009).Eelgrass metrics scaled with spacing between longlines due to both shading and dessication of eelgrass that coexisted, but stranded over the top of these lines in Humboldt Bay, California (55−65% cover with 3 m spacing reduced to <15% Researchers elsewhere found reduced eelgrass density, morphometrics, and photosynthesis due to reduced irradiance under suspended bags and oyster baskets, but also showed that the effect was mostly restricted to the area under the aquaculture structures themselves, and severity of effect de pended on the structure design and signature (Bulmer et al. 2012, Skinner et al. 2013, 2014).We suspect that the lack of an observed effect of longline culture at the landscape scale in our study was due to the very small area beneath these lines and inability to detect these changes at this broader scale.

Harvest method and temporal trajectory
Oyster aquaculture beds that were harvested with a mechanical dredge had significantly lower Z. marina than those harvested by hand or those reported to us as mixed harvest beds (sometimes harvested mechanically and sometimes by hand).Nonetheless, a mean of 99.9% of the model-predicted and 91.6% of the interpolation-predicted Z. marina were observed on mechanically harvested beds.The mean amount of Z. marina cover present in hand-harvested and mixed harvested beds was significantly higher and well over 100% of the predicted 40 Fig. 5. Relationship between the variability in Zostera marina present between years and the proportion of predicted Z. marina present for individual beds (top).When examined over time, individual beds exhibited 4 general patterns: (A) beds with chronically low Z. marina and low variation across years, (B) beds with the expected levels of Z. marina with low variation across years, (C) beds with highly variable levels of Z. marina across years, and (D) beds with high levels of Z. marina with high variation across years (lines represent trends in average predicted Z. marina/actual Z. marina, with red lines indicating 100%).Most beds cluster around B, with the expected amount of Z. marina cover remaining stable across years.All of the beds with lower than expected Z. marina present were mechanically harvest beds (x = 65%, n = 24) value, suggesting that culture actually enhanced the presence of eelgrass in these areas.These results on harvest impacts generally agree with previously collected data on wild shellfish harvest elsewhere (Peterson et al. 1987, Neckles et al. 2005) and in Willapa Bay, where mechanical harvest implements directly removed plants and caused more disturbance than hand-harvest or off-bottom longline oyster culture techniques (Tallis et al. 2009).Experimental use of a mechanical harvest dredge showed that both initial loss and the temporal pattern of recovery varied by site, with higher loss and slower recovery in areas with softer sediments (Tallis et al. 2009).Eelgrass seedling survival and growth was greater in mechanically harvested areas where more eelgrass had been removed and remaining shoots caused less shading (Wisehart et al. 2007).In our analyses, all of the oyster culture beds with < 65% of the mean expected amount of Z. marina cover (n = 24) were mechanically harvested beds and demonstrated a chronically low level of Z. marina present across all 3 years of study.While the majority of beds (all harvest types) had 65−145% of the model-predicted Z. marina and relatively low variability across years, the mechanically harvested beds with chronically lower amounts of eelgrass cover present were balanced by the presence of several mechanically harvested beds with very high (>145%) amounts of Z. marina cover present.Contrary to our expectation, the majority of beds with high year-to-year variability were hand-or mixed harvest beds and not mechanically harvested beds.Though our aquaculture spatial data was based on interviews conducted in a single year, we are confident that practices were fairly consistent over our 5 yr evaluation period, but this variability could be due to both aquaculture practices and other biotic and abiotic factors, such as eelgrass seedling recruitment success.Though presented in chronological order, the 3 years for which we collected data (2005,2006,2009) were not consecutive, and while trends were visible in the aerial photographs (Fig. 4), we did not examine actual harvest histories of individual beds by year, so it is difficult to precisely assess patterns of eelgrass loss and recovery.Experimental results from Willapa Bay suggest that recovery time for eelgrass after mechanical disturbance varied from 1 to 4 yr depending on initial impact and perhaps also on location within the estuary and eelgrass reproductive characteristics (i.e. the dispersal and success of seeds and seedlings versus vegetative propagation).Shading impacts from suspended culture and recovery times are less documented, but can also extend for more than a year (Skinner et al. 2014).
Further comparison and examination of data collected over multiple scales, including the intermediate scales of aquaculture beds and harvest cycles would be informative, but our results suggest that the majority of oyster aquaculture impacts are not persistent at the landscape scale.

Resilience, management implications, and research needs
Shellfish harvest and aquaculture have been important activities and supported local economies along the west coast of the US for well over 100 yr (MacKenzie et al. 1997, Dumbauld et al. 2011).Despite some negative results of disturbance to Z. marina from oyster aquaculture on relatively small spatial and short temporal scales, eelgrass is generally present and intermingles with oysters on all aquaculture beds at the tidal elevation where it is found in Willapa Bay, and our results suggest that current oyster aquaculture practices do not substantially reduce and may even enhance the presence of Z. marina at the estuarine landscape scale.Active culture of Crassostrea gigas currently occupies about 5% of the estuary (1764 ha) and 8% of the tideflat.With the exception of changes in practices such as switching from on-bottom culture to off-bottom culture in some locations, oyster culture disturbances have not changed materially for decades (Ruesink et al. 2006), so there is no reason eelgrass would necessarily be worse off now than in the recent past due to this disturbance.
Cultured C. gigas as habitat has likely historically replaced at least 3 other habitat types: monospecific eelgrass (Z.marina), habitat dominated by 2 species of burrowing thalassinid shrimp (Neotrypaea californiensis and Upogebia pugettensis), and native oyster Ostrea lurida habitat.Thom et al. (2003Thom et al. ( , 2014) ) provide evidence that eelgrass fluctuates with climate and environmental conditions, and there is compelling evidence that it has been expanding its distribution in Willapa Bay (Ruesink et al. 2010) and other estuaries along the open coast of the western US in the recent past, even though it is declining elsewhere in the world (Orth et al. 2006, Waycott et al. 2009) and in isolated locations on the US west coast (e.g.Hood Canal and San Juan Archipelago in Puget Sound; Gaeckle et al. 2007, Mumford 2007).Based on tidal elevation alone, Willapa Bay was estimated to contain 3139 ha of habitat suitable for Z. marina (0 to −1.2 m MLLW) in the 1850s, increasing to 4845 ha in the 1950s even as the bathymetry became shallower (Borde et al. 2003), but we found Z. marina well outside this tidal range.Here we estimate that Z. marina currently occupies about 6951 ha and Z. japonica covers at least an additional 3620 ha (B.R. Dumbauld unpubl.data).The expansion of Z. japonica since the mid 1990s and its potential facilitation of Z. marina (Bando 2006, Ruesink et al. 2012) offers one potential mechanism for the expanded distribution of Z. marina, but we have no data on shellfish culture as a third factor in this interaction.Burrowing shrimp also occupy a very large intertidal area in Willapa Bay (3060 ha, or 13.5% of the intertidal; B. R. Dumbauld et al. unpubl. data).Interestingly, the application of a pesticide to remove these shrimp (Dumbauld et al. 2006) is one aquaculture activity that could be responsible for enhanced eelgrass, because the shrimp negatively interact with eelgrass via sediment bioturbation and cause eelgrass seedlings to die (Dumbauld & Wyllie-Echeverria 2003, Berkenbusch et al. 2007).Shrimp have probably fluctuated in abundance and would have the potential to occupy a much larger area if shrimp control had not occurred.Finally, based on crude maps from the late 1800s, native oyster O. lurida beds occupied 6220 ha (17% of the low intertidal and shallow subtidal; Dumbauld et al. 2011, Zu Ermgassen et al. 2012) that now consists mostly of relatively undisturbed and dense Z. marina meadows.Thus, oyster habitat has always been present in Willapa Bay, although its current distribution and relief differ (e.g. C. gigas at higher tidal elevation and spread out instead of low intertidal reefs for O. lurida).Native oysters likely co-existed with eelgrass, stabilizing substrate and potentially creating pooled habitat and a mosaic at the landscape scale, just as C. gigas currently does when allowed to form reefs.
Eelgrass as habitat at the estuarine landscape appears to be resilient over both short and longer temporal periods and resistant to oyster aquaculture as a disturbance in this ecosystem.Resilience, originally defined by Holling (1973) as 'a measure of the ability to absorb changes of state variables, driving variables, and parameters and still persist' has been difficult for managers to address in part due to the lack of experimental tests and difficulty in defining measurable response thresholds (Hughes et al. 2005, Walker et al. 2006, Standish et al. 2014). Thom et al. (2012), however, defined resilience metrics and restoration actions that support resilient eelgrass habitat in US Pacific coast estuaries.These include sediment and water quality conditions at the landscape scale and shoot density at smaller spatial scales.Our research in Willapa Bay suggests that oyster aqua-culture as disturbance is generally within the scope of existing 'natural' disturbances to the system (e.g.winter storms), and eelgrass is inherently adapted to this scale of disturbance.Bivalve aquaculture has not been implicated in shifts to alternate states or reduced adaptive capacity of the larger ecological system.Typical thresholds that might be involved and have been linked to such catastrophic change and affect eelgrass at the landscape scale would likely be reached first with other human disturbances such as nutrients and eutrophication influencing the light environment or perhaps causing secondary shifts in macroalgal or herbivore abundance (Frederiksen et al. 2004, Orth et al. 2010b, Bostrom et al. 2014), or large-scale habitat removal and sediment alteration (Borde et al. 2003).Location and scale remain important management considerations, however, and our results from Willapa Bay may not represent other systems such as small inlets with stratified water columns and less routine physical disturbance by storms that could exhibit lower thresholds to aquaculture operations.
Although collective interest and emphasis on 'ecosystem-based management' and spatial planning at large scales has grown dramatically (Meffe et al. 2002, McLeod & Leslie 2009), in reality, management is driven by a wide variety of stakeholders and societal values (social, historical, political, moral, and esthetic, as well as economic), and progress towards these goals has been slow, especially in estuarine and marine systems (Arkema et al. 2006, Ruckelshaus et al. 2008, McLeod & Leslie 2009).Since seagrasses are widely known to provide valuable nursery habitat for fish and invertebrates in estuaries worldwide (Bostrom et al. 2006, Gillanders 2007, Heck et al. 2008), and in particular for juvenile salmon on the west coast of the US (Aitkin 1998, Fresh 2006), most management policy is still based on a general 'no net loss' rule for all estuarine wetlands.Broader management policies have been codified in the US National Ocean Policy (US Nation Ocean Council 2013), and marine spatial planning efforts and some nascent integrated ecosystem assessment efforts are underway (Levin et al. 2009, Foley et al. 2010, Samhouri et al. 2014), but most permitting actions, including those for aquaculture, are still conducted at the scale of local and individual actions.The US ACOE issued nationwide permit 48 in 2007, which covered existing shellfish aquaculture operations, and continues to work with the shellfish aquaculture industry and other federal and state agencies to renew these permits and implement this permitting and reporting program, which includes permit-ting actions for new aquaculture activities and operations.Recent listings of several stocks of salmon and green sturgeon under the US Endangered Species Act (Adams et al. 2002, Good et al. 2007) and revisions in the Magnuson Stevens Act designed to protect essential fish habitat have required both aquaculture industry and US ACOE to consult with other agencies, including NOAA, on these actions.Existing aquaculture activities were generally allowed to continue, but 'pre-construction notification' was required if mechanical harvest, tilling, or harrowing was conducted in areas with eelgrass present and individual permits required if culture was to be allowed in a new project area with > 0.2 ha area of eelgrass present.A very precautionary approach has been taken on these new permits, and interactions with Z. marina mostly prohibited at very small scales (e.g.including 5 m buffer zones from edges of eelgrass meadows, which are defined based on data collected at the time of permit issuance and thus presence in a single year in the estuary where they occur).Mitigation is required for all loss and little to no consideration made for location and scale or for the value of other estuarine habitats (oysters or open mudflat).Our data suggest that at least in Willapa Bay, a broader ecosystem perspective would be useful and could be less restrictive with regards to protecting eelgrass per se, though further research should be conducted to determine the functional value of eelgrass and shellfish as habitat for the resources being protected at this larger scale.
Use of eelgrass and shellfish aquaculture as habitat and the behavioral response of commercially important species (salmon, crab) have been studied across habitat types in Willapa Bay and a few other US west coast estuaries (Holsman et al. 2006, Hosack et al. 2006, Semmens 2008, Pinnix et al. 2013, Dumbauld et al. in press); however, these studies have mostly examined organism presence and density in a given habitat and not broad-scale spatial pattern or functional roles of these habitats.The influence of habitat configuration and human changes to that configuration on organism abundance and behavior at broad spatial scales (relative to the organism being studied) have been widely examined in terrestrial systems (Debinski & Holt 2000, Kareiva et al. 2007, Lindenmayer et al. 2008), but only recently considered for marine habitats such as eelgrass meadows (Bostrom et al. 2006, Connolly & Hindell 2006, Hinchey et al. 2008) and oyster reefs (Harwell 2004, Grabowski et al. 2005).It could be that some habitats are more important than others at a broader landscape scale (e.g.eelgrass as protective cover near channel edges for juvenile salmon; Semmens 2008) and that the configuration of both oysters and eelgrass as habitat is also important because it provides food for larger, more mobile organisms at that scale (e.g. for juvenile salmon, English sole, or shorebirds and waterfowl) and protective cover and food for others (e.g. for juvenile crab; Fernandez et al. 1993, Tait & Hovel 2012).Though less is known about habitat function for these larger, more mobile organisms, they can use structure for protection from even larger predators (juvenile salmon in eelgrass and 0+ Dungeness crab in oyster), but still rely on other habitats for foraging (e.g.1+ juvenile Dungeness crab in open unstructured habitat; Holsman et al. 2006).Given the presence of mixed habitats (i.e.eelgrass within aquaculture beds), new research should address landscape-level features such as patch size, connectivity, and the population response of these organisms to the estuarine habitat mosaic.Fragmentation, patchy seagrass, and increased habitat edges may actually enhance diversity and increase the density of some bottom-feeding invertebrates such as decapod crustaceans and fish, whereas larger seagrass meadows may harbor higher numbers of smaller cryptic species (Salita et al. 2003, Tanner 2005, Selgrath et al. 2007, Horinouchi et al. 2009).Focus on these broader landscape patterns instead of individual actions and use at small scales would also be an area where best management plans could be designed and implemented, since the shellfish industry would likely be supportive of maintaining habitat corridors (e.g.eelgrass along channel edges) and timing windows (e.g.limited harvest operations in some areas during fish spawning or bird migration periods), should this prove advantageous to resource management at this scale.

Fig. 1 .
Fig. 1.Willapa Bay estuary.Note large intertidal area of the estuary (62%, light gray).Insets -top right: location in Washington State, USA; middle and lower right: shellfish aquaculture beds (dark grey) and eelgrass Zostera marina (black) distribution in 2005 in a northeastern portion of the estuary

Fig. 2 .
Fig. 2. (a) Frequency distribution of the area of intertidal zone, Zostera marina, and oyster aquaculture by tidal elevation and (b) distance to the estuary mouth in Willapa Bay.MLLW: mean lower low water (small bars for eelgrass represent SE)

Fig. 3 .
Fig. 3. (a) Relationship between predicted and observed eelgrass Zostera marina cover in all aquaculture beds for 2005 and (b) relationship between interpolated and observed Z. marina in all aquaculture beds for 2005.Model fit lines and confidence intervals overlap the 1:1 line, which indicates they had the same amount of Z. marina cover present as predicted; however, a group of beds fell well below and outside the 95% confidence and prediction intervals

Fig. 4 .
Fig. 4. Aerial infrared photographs of tideflats in northern Willapa Bay from 2005, 2006, and 2009, with oyster aquaculture beds outlined in black.Note changes over time, with vegetation (red) evident on some beds (A) in 2005, absent on others (B), and changing over time in 2006 and 2009.Vegetation also changed in areas outside aquaculture beds (C) fold cross validation RMSE = 0.22, 0.23, 0.17, respectively).Models ex ploring individual predictors suggest that distance to estuary mouth and tidal elevation explained most of the variation in Z. marina cover.The best-performing models only describe approximately half of the variation in Z. marina cover in each year(2005: 44.5%, 2006: 41.2%, and  2009: 53.5%;Table2;TableS1in the Supplement).

Table 1 .
Area statistics for Willapa Bay estuary.Active intertidal oyster aquaculture excludes clam beds, seed racks, subtidal beds, and partially active beds.Vegetation includes all vegetation signal between −1.5 m and + 3.2 m mean lower low water (MLLW) and includes some Spartina alterniflora, dense Zostera japonica, and Z. marina in pools at higher elevation, calculated using the mean of the summed probabilities from2005, 2006, and 2009.Percentages of sites occurring at, below, or above 1.2 m MLLW are given for 4238 ground survey sites with 3 different densities of Z. marina cover and for mapped S.
a File received from the Olympic Natural Resources Center, Forks, WA

Table 2
. Relative contribution of predictor variables to the generalized additive model (GAM) model for Zostera marina in 2005 (see Table S1 in the Supplement for 2006 and 2009 results).Contributions of individual predictors are shown for a model with smoothing for each predictor.The final model incorporates a tensor smoothe including 3 predictors and individual smoothes for 2 predictors.Deviance explained is given in parentheses.edf = estimated degrees of freedom, χ 2 value reflects the relative importance of each predictor in the model.The p-values are provided but are only useful for identifying predictors that do not contribute to the model

Table 3
. Total estimated area and proportion calculations for Zostera marina and oyster aquaculture in Willapa Bay in all 3 years using both model and interpolation predictions

Table 4 .
Mixed effects model comparison using maximum log-likelihood to determine the amount that harvest method (hand, dredged, mixed) and oyster culture bed type (seed, fattening, direct, longlines) contribute to explaining the variation in the actual to model-predicted and interpolation-predicted amounts of Zostera marina cover on each bed across the 3 observation years.Bed ID and Year are included as random effects.The model including only harvest method has a lower Bayesian information criterion (BIC) score in both cases.na: not applicable