ABSTRACT Title of Document: POPULATION DYNAMICS OF EASTERN OYSTERS (CRASSOSTREA VIRGINICA) IN THE CHOPTANK RIVER COMPLEX, MARYLAND, DURING 1989 - 2015 Matthew Daniel Damiano, Master of Science, 2017 Directed By: Professor Michael J. Wilberg Marine Estuarine Environmental Sciences The eastern oyster (Crassostrea virginica) fishery in the Choptank River Complex (CRC) supports a large fraction of Maryland’s harvest. The CRC is also host to some of the largest oyster restoration projects in the world. Yet the relative effects of harvest and restoration on the population dynamics of oysters in the CRC have not been assessed. We developed stage-based population models for each region of the CRC in AD Model Builder using dredge survey and harvest data provided by Maryland Department of Natural Resources from 1989 to 2015. Natural mortality was low during 2004-2015 potentially due to increased resistance to the disease, Dermo. Recruitment was greatest in the late 1990s, 2010, and 2012, which caused an increase in abundance. These models will serve as the foundation of a simulation model that will be used to help fishery stakeholders evaluate management and restoration options in the CRC. POPULATION DYNAMICS OF EASTERN OYSTERS (CRASSOSTREA VIRGINICA) IN THE CHOPTANK RIVER COMPLEX, MARYLAND, DURING 1989 - 2015 By Matthew Daniel Damiano Thesis submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment of the requirements for the degree of Master of Science 2017 Advisory Committee: Professor Michael J. Wilberg, Chair Assistant Research Professor Dong Liang Associate Professor Elizabeth North © Copyright by Matthew Daniel Damiano 2017 ii Acknowledgements I would like to acknowledge and thank the National Science Foundation for funding my research, and Maryland Department of Natural Resources, the National Oceanic and Atmospheric Administration, the Oyster Recovery Partnership, and University of Maryland Center for Environmental Science Paynter Lab for sharing their data. Special thanks to the OysterFutures work group for their insight and assistance throughout the duration of my work. Without these groups, none of this would have been possible. I would like to thank my family for their unflagging love and support. Amber Walker has been instrumental in helping me navigate the challenges of balancing work and life. She routinely inspired me to work hard, and be strong in the face of any obstacle. My mother, Anne Tileston, father, Robert Damiano, and brother, William Damiano have supported my passion for marine science from that first spark of interest in the tide pools of Maine, and I would not have made it this far without them. I would also like to acknowledge Mike Wilberg, Elizabeth North, Genny Nesslage, and Tom Miller for their strong support and professional guidance. Additionally, I thank the Wilberg Lab for their feedback throughout the last two years. To my friends, Brian Gallagher, Zoraida Paola Perez-Delgado, Kathryn Doering, and Suzan Shahrestani, your company, collegiality, and fierce intellect have sustained me throughout this journey. iii Table of Contents Acknowledgements .......................................................................................................................... ii Table of Contents ............................................................................................................................ iii Table of Tables ................................................................................................................................ v Table of Figures .............................................................................................................................. vi Introduction .................................................................................................................................... 1 Methods ........................................................................................................................................... 5 Overview ...................................................................................................................................... 5 Data ............................................................................................................................................. 7 1. MD DNR buy tickets ....................................................................................................... 7 2. Maryland fall dredge survey ............................................................................................ 9 3. Sanctuary monitoring ..................................................................................................... 10 4. Patent tong population estimates.................................................................................... 10 5. Spat on shell planting ..................................................................................................... 11 6. Habitat restoration .......................................................................................................... 11 7. SONAR bottom mapping ............................................................................................... 12 Assessment model ...................................................................................................................... 12 1. Population model ........................................................................................................... 12 2. Observation model ......................................................................................................... 15 3. Model fitting .................................................................................................................. 16 Sensitivity analyses .................................................................................................................... 20 Correlation tests..................................................................................................................... 20 Results ........................................................................................................................................... 21 Model fits ................................................................................................................................... 21 Mortality .................................................................................................................................... 26 Recruitment ................................................................................................................................ 27 Abundance ................................................................................................................................. 28 Habitat ....................................................................................................................................... 29 Exploitation ................................................................................................................................ 30 Parameter estimates................................................................................................................... 32 Sensitivity analyses .................................................................................................................... 33 Correlation tests ........................................................................................................................ 34 iv Discussion ..................................................................................................................................... 36 Appendices .................................................................................................................................... 43 Appendix 1. Table of spat on shell planting by the Oyster Recovery Partnership (ORP) by region. Regions are identified using abbreviations: Broad Creek (BC), Lower Choptank River (LC), Middle Choptank River (MC), Upper Choptank River (UC), Tred Avon River (TA), Little Choptank Open (LCO), Little Choptank Closed (LCC), Harris Creek Open (HCO), and Harris Creek Closed (HCC). ...................................................................................................... 43 References ..................................................................................................................................... 45 v Table of Tables Table 1. Definition of model variables and parameters. ................................................................ 18 Table 2. Estimated values of parameters that were assumed constant in the model. G represents the probability of an oyster growing from the small to market stage (per year), h the rate of habitat loss (% per year), d the instantaneous rate of box disarticulation (per year), and q the catchability coefficient. The subscripts for d and q indicate the stage: spat (sp), small (sm), market (mk), small box (smb), and market box (mkb). Regions are identified using the following abbreviations: Broad Creek (BC), Lower Choptank River (LC), Middle Choptank River (MC), Upper Choptank River (UC), Tred Avon River (TA), Little Choptank Open (LCO), Little Choptank Closed (LCC), Harris Creek Open (HCO), and Harris Creek Closed (HCC). The subscripts for d and q indicate the stage: spat (sp), small (sm), market (mk), small box (smb), and market box (mkb). .......................................................................................................................... 33 Table 3. Sensitivity analyses results for the Broad Creek estimation model. Values indicate the relative differences (%) from the parameter values obtained from the base model. Parameter values are from left to right: growth (G), average annual rate of natural mortality (M), average annual recruitment (r), and average annual adult abundance (N), and rate of habitat loss (h). Please, refer to the methods section for scenario descriptions. ...................................................... 34 Table 4. Pearson correlation coefficients (above the diagonal) and p-values (shaded region below the diagonal) for the natural mortality rate among regions. Regions are identified using abbreviations: Broad Creek (BC), Lower Choptank River (LC), Middle Choptank River (MC), Upper Choptank River (UC), Tred Avon River (TA), Little Choptank Open (LCO), Little Choptank Closed (LCC), Harris Creek Open (HCO), and Harris Creek Closed (HCC). .............. 35 Table 5. Pearson correlation coefficients (above the diagonal) and p-values (shaded region below the diagonal) for recruitment among regions. Regions are identified using abbreviations: Broad Creek (BC), Lower Choptank River (LC), Middle Choptank River (MC), Upper Choptank River (UC), Tred Avon River (TA), Little Choptank Open (LCO), Little Choptank Closed (LCC), Harris Creek Open (HCO), and Harris Creek Closed (HCC). ....................................................... 36 vi Table of Figures Figure 1. Map of the CRC by NOAA code: 053 (Little Choptank), 137 (Lower Choptank), 237 (Middle Choptank), 337 (Upper Choptank), 437 (Harris Creek), 537 (Broad Creek), 637 (Tred Avon). .............................................................................................................................................. 7 Figure 2. Map of regions of the CRC that are open (tan) and closed (green) to harvest. ................ 7 Figure 3. Diagram of the fate of oysters in each stage of the population model. ........................... 14 Figure 4. Estimated (line) and observed (points) log relative density of spat in number of oysters per bushel of cultch from the MD DNR fall dredge survey in all regions of the Choptank River Complex during 1989-2015. .......................................................................................................... 21 Figure 5. Estimated (line) log relative density of small oysters in number of oysters per bushel of cultch fit to observed density (points) from the MD DNR fall dredge survey in all regions of the Choptank River Complex from 1989-2015. .................................................................................. 22 Figure 6. Estimated (line) log relative density of market oysters in number of oysters per bushel of cultch fit to observed density (points) from the MD DNR fall dredge survey in all regions of the Choptank River Complex during 1989-2015. .......................................................................... 22 Figure 7. Estimated (line) log relative density of small boxes in number of boxes per bushel of cultch fit to observed density (points) from the MD DNR fall dredge survey in all regions of the Choptank River Complex during 1989-2015. ................................................................................ 23 Figure 8. Estimated (line) log relative density of market boxes in number of boxes per bushel of cultch fit to observed density (points) from the MD DNR fall dredge survey in all regions of the Choptank River Complex during 1989-2015. ................................................................................ 23 Figure 9. Estimated (line) and observed (points) of log catch per unit effort (CPUE) for the hand tong fishery at the beginning of the season (October) from the MD DNR buy ticket data for all regions of the Choptank River Complex during 1989-2015. ......................................................... 24 Figure 10. Estimated (line) and observed (points) of log catch per unit effort (CPUE) for the hand tong fishery at the end of the season (March) from the MD DNR buy ticket data for all regions of the Choptank River Complex during 1989-2015. .......................................................................... 25 Figure 11. Estimated (line) and observed (points) of log catch per unit effort (CPUE) for the power dredge fishery at the beginning of the season (November) from the MD DNR buy ticket data for all regions of the Choptank River Complex during 1989-2015 ........................................ 25 Figure 12. Estimated (line) and observed (points) of log catch per unit effort (CPUE) for the power dredge fishery at the end of the season (March) from the MD DNR buy ticket data for all regions of the Choptank River Complex during 1989-2015. ......................................................... 26 Figure 13. Estimated natural mortality (% per year) for small and market adult oysters and 95% confidence intervals (shaded) in each region of the Choptank River Complex during 1989 - 2015. .............................................................................................................................................. 27 Figure 14. Estimated recruitment with 95% confidence intervals (shaded) in each region of the Choptank River Complex during 1989 – 2015. ............................................................................. 28 Figure 15. Estimated adult (small and market) eastern oyster abundance with 95% confidence intervals (shaded) in each region of the Choptank River Complex during 1989 – 2015. .............. 29 Figure 16. Estimated change in the area of hard bottom oyster habitat and 95% confidence intervals for each region of the Choptank River Complex during 1989 - 2015. ............................ 30 vii Figure 17. Estimated exploitation rate (% per year) of market-sized eastern oysters with 95% confidence intervals (shaded) for all regions of the Choptank River Complex during 1989 - 2015. .............................................................................................................................................. 31 Figure 18. Total annual mortality from fishing and natural mortality. Percent of oysters lost to exploitation is white under the shaded total (% of oysters lost to all sources of mortality). ......... 32 1 Introduction Oysters are bivalve mollusks, and many species are ecosystem engineers that form reef habitat through aggregation of conspecifics (Coen et al., 1999). These reefs provide biogenic habitat for conspecifics and numerous marine and estuarine species of reef-associated fishes and invertebrates (Lenihan et al., 1999; Beck et al., 2011). Oysters’ reliance on oyster shell for habitat make them particularly vulnerable to the effects of fishing (Wilberg et al., 2013), and it is estimated that there has been an 85% reduction in oyster reef ecosystems worldwide (Beck et al., 2011). Eastern oysters (Crassostrea virginica) are a commercially and ecologically important species that inhabit coastal waters and estuaries along the Atlantic and Gulf Coasts of North America. They provide ecosystem services such as improving water quality and clarity through seston reduction, providing complex hard-bottom habitat, and promoting biodiversity (Newell, 1988; Coen et al., 2007). Of these benefits, the eastern oyster’s ability to reduce solids and phytoplankton in the water column through suspension feeding has been observed in situ in shallow, mesohaline waters in Chesapeake Bay; although some of the other potential benefits are not supported by the same degree of empirical evidence (Coen et al., 2007). Eastern oyster reefs promote biodiversity in several ways: they provide refuge to fishes and invertebrates, act as a coupling between benthic and pelagic systems, and bolster fishery production (Lenihan & Peterson, 1998). The Maryland fishery for eastern oysters in Chesapeake Bay has experienced a large decline due to overharvest, loss of habitat, and several outbreaks of disease (Rothschild et al. 1994, Wilberg et al. 2011). The fishery reached peak harvest in 1884 at 615,000 metric tons (Rothschild et al., 1994). Annual harvest has since dropped to approximately 3% of the fishery’s 2 peak (Tarnowski, 2016). More efficient gear types such as dredge vessels and hydraulic patent tongs have been introduced over time that allowed the Maryland fleet to expand into deeper waters and remove eastern oysters at a higher rate (Rothschild et al., 1994). Although eastern oyster management has restricted harvest to the use of particular gear types in certain regions, the use of power dredging, a more efficient dredge gear (Chai et al., 1992) expanded into upper Chesapeake Bay in 2003. In addition, an overall increase in fishing activity drove the rate of exploitation from 5% in 2003 to an average of 25% during 2004 – 2008 (Wilberg et al., 2011). This is above the estimated range of the exploitation rate that would produce maximum sustainable yield (5 – 10%) for Chesapeake Bay, Maryland (Wilberg et al., 2013). Abundance of eastern oysters in Maryland was estimated at 0.3% of the abundance before the onset of commercial fishing (Wilberg et al., 2011). The recent decline in oyster populations has also been attributed to increased natural mortality from the introduction of two diseases into Chesapeake Bay (Mann & Powell, 2007).The diseases MSX (Haplosporidium nelsoni) and Dermo (Perkinsus marinus) both caused extensive epizootic mortality events in eastern oysters in Chesapeake Bay during the 1960s in high salinity waters (Burreson et al., 2000). The spread of MSX is attributed to the introduction of infected Japanese oysters (Crassostrea gigas), while the origins of Dermo are more nebulous (Mann & Powell, 2007), and Dermo may have been associated with eastern oysters in Chesapeake Bay before its discovery (Burreson & Ragon Calvo, 1996). Although both diseases tended to be found in eastern oysters in areas with higher salinity, during years of drought, 1986 - 1987 and 1990 - 1992, Dermo expanded its range into the mesohaline parts of Chesapeake Bay, Maryland (Cook et al., 1998). Estimates of annual natural mortality during 1986 – 1987 were approximately 60%, and 25% during 1990 – 1992 in Maryland, albeit with substantial 3 uncertainty (Wilberg et al., 2011). Dermo remains prevalent among eastern oyster populations in Maryland (Tarnowski, 2016), but large mortality events have not occurred since 2002 – 2003 when the average natural mortality rate was close to 50% (Wilberg et al., 2011). Management and restoration efforts have been implemented to counter the effects of overharvest and habitat loss. In the late 1800s, the Maryland Oyster Commission and the Maryland Oyster Police were formed to improve and protect the fishery. In the early 20th century, a minimum size limit of three inches was imposed on the fishery as well as a shell tax to provide shell for replenishment activities (Kennedy & Breisch, 1983). During 1961 – 2016, Maryland established 29 oyster sanctuaries, with the three largest sanctuaries established in 2010 in the Choptank River Complex (Maryland Dept. of Natural Resources, 2016). The Choptank River Complex sanctuaries were designed as part of an effort to protect 24% of the remaining oyster resource in Maryland. Additional restoration efforts have largely involved deploying young-of-the-year oysters (spat) on shell in depleted oyster habitat (Mann & Powell, 2007), although artificial reefs have been constructed in the three largest oyster sanctuaries since 2010 (Westby et al., 2016). Implementing restoration and management strategies to counter the decline has been a contentious issue between natural resource managers and oyster harvesters (Bocking, 2011). In response, the OysterFutures facilitated workgroup was formed in summer 2015 with the goal of developing consensus solutions for oyster restoration and management in the Choptank River Complex (CRC); the CRC is the site of the largest oyster restoration effort in the United States (Westby et al. 2016) and produces 28% of the annual eastern oyster harvest in Maryland on average (Tarnowski 2016). OysterFutures is a collaborative, facilitated process in which scientists work with stakeholders to determine fishery management objectives, develop 4 strategies to achieve those objectives, and present them to the management agency. OysterFutures follows a similar format as the FishSmart process (Ihde et al., 2010). The OysterFutures work group is composed of stakeholders in the Maryland oyster fishery including watermen, science educators, non-governmental conservation organizations, aquaculturists, agency representatives, and oyster buyers who are advised by a panel of scientists and facilitators from the University of Maryland Center for Environmental Science, Florida State University, and the Virginia Institute of Marine Science. As part of the process, the work group is collaboratively building a model to forecast the potential effects of alternative restoration and management options. The workgroup will conclude by making management recommendations to Maryland Department of Natural Resources (MD DNR). Improved understanding of the population dynamics of eastern oysters in the CRC is important for several reasons. The CRC represents a substantial component of the Maryland fishery at 28% of annual yield (Tarnowski, 2016). However, management of oysters in Maryland does not currently use abundance or fishing mortality rate estimates in determining management measures. A framework in which eastern oyster can be assessed at the tributary- level may provide valuable information for management of oysters in this region. Furthermore, efforts to conduct stock assessments for all of Maryland will benefit from the development of methods to estimate abundance and mortality rates in the CRC. A population dynamics approach can also be useful for evaluating the short-term efficacy of sanctuaries in mesohaline tributaries. Lastly, models include general fishery processes such as fishing, natural mortality, recruitment, and growth, and could potentially be applicable to other commercially exploited bivalves outside of the CRC. 5 Although assessments have been performed to evaluate the effects of fishing and disease on oysters in the Maryland portion of Chesapeake Bay (Rothschild et al. 1994, Wilberg at al. 2011), the degree to which harvest, restoration, and natural mortality (including disease) have affected the abundance of oysters in the CRC has not yet been assessed. Nor are current population estimates available for the CRC. Therefore, this thesis is driven by two primary objectives: 1) to estimate abundance and rates of natural mortality and exploitation in all regions of the CRC, and 2) to compare estimates of natural mortality rates and recruitment among regions. To achieve these objectives, we implemented stock assessment models for each of the seven regions of the CRC to estimate how abundance and rates of natural mortality and exploitation have changed over time. Models were fitted to data from the MD DNR fall dredge survey, fishery catch per unit effort (CPUE) from the MD DNR buy ticket program, and estimated abundance derived from monitoring activities in Broad Creek, Harris Creek, and the Little Choptank River. Methods Overview We built statistical, stage-structured models for oysters in the CRC similar to the one developed by Wilberg et al. (2011). Models were developed in AD Model Builder (Fournier et al., 2012) to estimate abundance by stage, natural mortality rates, amount of available habitat, and exploitation rates. The models were fitted to indices of density from the MD DNR fall dredge survey and fishery catch per unit effort (CPUE) from dealer reported buy tickets during 1988 – 2016. Models were developed for nine regions of the CRC. The model year began on October 1st, corresponding to the beginning of the fishing season. Years in the model refer to the 6 year of the end of the fishing season; for example, the model year for October 1, 1988 through September 30, 1989 were labeled as 1989. Study Area The CRC is an estuary and tributary containing several smaller tributaries located in upper Chesapeake Bay, Maryland. The CRC is approximately 2000 km2 with a salinity range of 5 to 20 depending on the location and time of year. The CRC was divided into seven regions by harvest reporting code (the lowest level of reporting for harvest and effort): Little Choptank River (053), Lower Choptank River (137), Middle Choptank River (237), Upper Choptank River (337), Harris Creek (437), Broad Creek (537), and Tred Avon River (637) (Figure 1). Several regions contain areas closed to harvest for human health or as oyster sanctuaries. For these two of these regions, Harris Creek and the Little Choptank River, areas open and closed to harvest were modeled separately to improve the accuracy of estimation (Figure 2; Pincin & Wilberg, 2012). 7 Figure 1. Map of the CRC by NOAA code: 053 (Little Choptank), 137 (Lower Choptank), 237 (Middle Choptank), 337 (Upper Choptank), 437 (Harris Creek), 537 (Broad Creek), 637 (Tred Avon). Figure 2. Map of regions of the CRC that are open (tan) and closed (green) to harvest. Data The data included amount of harvest from dealer reported buy tickets and counts of oysters from the fall dredge survey provided by MD DNR, patent tong monitoring data, oyster spat on shell planting data from the Oyster Recovery Partnership (ORP), and habitat restoration data and estimates of available oyster habitat from SONAR bottom mapping provided by NOAA. 1. MD DNR buy tickets Oyster buyers are required to submit reports called buy tickets to MD DNR that include the number of bushels of oysters that were purchased, hours spent harvesting, location of harvest, and the number of individuals on the vessel (Tarnowski, 2016). Raw data from the buy tickets were available for 1989 – 2015 and were used to calculate total harvest by region in number bushels. The location of harvest is reported by NOAA code. 8 Bushels of oysters harvested were converted to numbers of oysters harvested. A Maryland bushel for oyster harvest is slightly larger than the US agricultural bushel at 45.9 liters. We used an average number of oyster per Maryland bushel of 258 oysters at or above the legal size limit (76 mm shell height) based on estimates from oyster buyers and harvesters in. Up to 5% of a bushel may be composed of sub-legal oysters (Tarnowski, 2016). We extrapolated the number of smalls per bushel, 360 smalls per bushel, from estimates of the number of oysters per bushel at different sizes) using an average size for an undersized oyster of 64 mm. We assumed that 5% of the harvest by volume was sub-legal based on discussions with harvesters and buyers. We corrected the harvest estimates for underreporting. Conversations with MD DNR scientists and commercial watermen suggested that the harvest data were underestimates of the true amount harvested. Prior to 2010 harvest data were only provided by oyster buyers, but in 2010 harvesters were required to report daily harvest and effort. We assumed a rate of 24% underreporting for years prior to 2010 based on discussions with MD DNR (Frank Marenghi, MD DNR, personal comm.). Percent rates of reporting from harvester reports were provided by MD DNR for 2010 through 2015. On average, 86% of expected harvest reports were submitted during 2010 - 2015; we assumed that non-reporting was random and used 14% as an estimated underreporting rate for 2010-2015 (Frank Marenghi, MD DNR, personal comm.). We corrected the time series of harvest in bushels using these underreporting rates. We calculated catch per unit effort (CPUE; bushels per man hour) for the hand tong and power dredge fisheries at the beginning and end of the season to provide the models with information on the amount of depletion during the fishing season. Hand tong and power dredge gears accounted for 91% of landings in the CRC during 1989 – 2015. Power dredging is more efficient at removing oysters than other gears (Chai et al. 1992), so we calculated separate CPUE 9 indices for the two gears. CPUE was calculated as the geometric mean number of bushels harvested per person hour for each year by gear type from the first month of harvest (October for hand tong, November for power dredge), and the last month of harvest (March). Power dredging was expanded in the CRC in 2003 (Tarnowski, 2016), and became the dominant means of harvest in the Lower Choptank, and Harris Creek Open regions. CPUE data were not always available for each month and gear either because information was missing on the buy ticket to calculate effort or no fishing occurred during that period for a gear. For years with CPUE data available for the beginning or end of the season, only one index was calculated. 2. Maryland fall dredge survey Fall dredge survey data were available from 1989 – 2015. The fall dredge survey has been conducted annually from October through December using a motorized vessel that deployed and towed an 81-cm-wide oyster dredge over 259 natural oyster bars, including public fishery bars and bars located within sanctuaries. A 1/2 Maryland bushel sub-sample was taken from the dredged material (cultch) after each tow. The half-bushel subsample was sorted by stage and counted. The stages included young-of-the-year oysters (spat), small adult oysters (smalls; <76 mm shell height), and market-sized adult oysters (markets; >76 mm shell height). Additionally, the articulated shells (valves still attached at the hinge ligament) of dead oysters (boxes) were recorded for the same size categories as live oysters (smalls, markets; Jordan et al., 2002; Tarnowski, 2016). Spat boxes were also recorded but were very rare due to their fragility. All counts were then recorded as number per Maryland bushel. We used the Fall Dredge Survey data to develop indices of relative density (Wilberg et al. 2011) for each stage. We fitted generalized linear mixed effect models (GLMM) to the Fall Dredge Survey count data (Maunder & Punt, 2004; Bolker et al., 2009). The GLMMs included 10 year as a fixed effect (βy), and site (i.e., oyster reef; γbar) as a random effect (Eq. 1) to correct for changes in sites sampled over time and potential differences in catchability among sites (Wilberg et al. 2011), 𝑙𝑜𝑔𝑒(𝐼𝑦,𝑠) = 𝛽𝑦 + 𝛾𝑏𝑎𝑟 + ԑ. (1) The model included a log link function and a negative binomial distribution to account for the over-dispersion of counts in Fall Dredge Survey (Linden & Mantyniemi, 2011). The model was applied separately to data for each region to develop indices of relative density by stage. 3. Sanctuary monitoring Patent tong surveys were conducted in Harris Creek to monitor progress of restoration efforts in that region. The Harris Creek Closed model used patent tong monitoring data collected during 2012 – 2015. These data were available at the bar level in number of oysters caught in a patent tong grab for both restored and unrestored habitat and on five oyster bars monitored annually since 2012. Density (number per m2) values by stage of live oysters were calculated by dividing the number of oysters caught by the area of the patent tong. To obtain an estimate of abundance for live oysters (spat, smalls, and markets) in 2015, live densities by stage were averaged for restored and unrestored habitat and multiplied by an estimate of total amount of each habitat. 4. Patent tong population estimates Data were available for Harris Creek, Broad Creek and the Little Choptank from patent tong surveys conducted to validate the side-scan SONAR bottom maps in 2010 and provide estimates of oyster population size and density. Data were collected during summer months using a stratified random sampling approach with strata classified by type and amount of available habitat. 30 to 314 points per stratum were surveyed. Estimates of abundance in each habitat classification were calculated as the product of mean density and the area of each habitat 11 classification. The total abundance of oysters in a region was the sum of the abundance estimate in each habitat classification. Bootstrapping was used to estimate the precision of the estimates. The estimated number of live oysters (spat, smalls, and markets) were 10 million for the Little Choptank Closed region, and 23 million for the Harris Creek Closed region (K. Paynter, University of Maryland Center for Environmental Science, unpublished data). 5. Spat on shell planting ORP has planted hatchery-reared spat-on-shell as part of restoration activities in the CRC since 1998. ORP provided data on the number of spat planted by region (Appendix 1). Approximately one-month old spat that were 1-2 mm shell height were planted (Paynter et al., 2010). The number of eastern oysters planted varied by region and year. During 1998 – 2008, the Middle Choptank was seeded with 15 million spat, and the Upper Choptank with 460 million (ORP, unpublished data). During 2011 – 2015, Harris Creek has been planted with 2 billion spat, the Little Choptank with 130 million, and the Tred Avon with 10 million (Westby et al., 2016). Survival of spat from planting to October in the model was assumed to be 15% based on estimated survival from dive surveys in Harris Creek (Paynter et al., 2014) 6. Habitat restoration Major restoration activities began in the CRC in 2010 in which oyster reefs were restored by planting hatchery-reared spat-on-shell and rebuilding hard bottom substrate using stone, clam shell, and oyster shell. As of the last year used in the models (2015), 0.55 km2 of oyster reefs had been constructed from rock and bivalve shell and seeded with spat in Harris Creek, 0.52 km2 in the Little Choptank River, and 0.06 km2 in the Tred Avon (Westby et al., 2016). Estimates of repletion activities by MD DNR using fossil shell during 1984 – 1999 were not included in the models because these habitat plantings comprised approximately 0.02% of the estimated 90 km2 12 of hard bottom habitat in the CRC as of 2001 (Smith et al., 2005), 0.03% of the estimated 51 km2 of habitat in the CRC as of 2010 (Allen et al., 2013). Shell from these activities was rapidly lost to sedimentation within about five years after planting (Smith et al., 2005). 7. SONAR bottom mapping The area of available hard bottom habitat for oysters in the CRC was estimated by harvest reporting region in 2010. Side-scan SONAR surveys were performed by the Maryland Geologic Survey and NOAA. NOAA created GIS polygons of oyster habitat by combining the SONAR data with patent tong survey data, ponar grabs, and acoustic classifications (Allen et al., 2013). Assessment model 1. Population model The models were stage-based using the five stages from the fall dredge survey: spat, smalls, markets, small boxes, and market boxes (Figure 3). The model year begins on October 1, at the beginning of the hang tong season, and coincides with the timing of the fall dredge survey. The processes being modeled include recruitment to the fishery (natural and planted), growth from small to market, natural mortality (that includes disease) of spat, smalls and markets, the effect of fishing on small and market oysters (fishing mortality), changes to habitat over time, and the disarticulation of small and market boxes. Model variables and parameters are described in Table 1. The initial abundance of spat, smalls, and markets in the first year were estimated parameters. The model estimated the spat abundance each year as the sum of the estimated number of naturally produced spat and the number of planted spat (Eq. 2), 𝑁𝑦+1,0 = 𝑒 ?̅?+?̌?𝑦 + 𝑁𝑦,𝑎𝑆𝑦,𝑎. (2) 13 The number of spat planted was multiplied by a survival rate of 15% based on densities of spat one to two months after planting relative to the initial planting density (K. Paynter, UMCES, unpublished data). The number of smalls each year was estimated by calculating the number of spat that survive natural mortality to become smalls and adding those to the smalls already in the population that survived natural mortality and harvest, but did not grow to be markets (Eq. 3), 𝑁𝑦+1,1 = (𝑁𝑦,1 − 𝐶𝑦,1)(1 − 𝐺)𝑒 −𝑀𝑦 + 𝑁𝑦,0𝑒 −𝑀𝑠𝑝. (3) Harvest was modeled as a pulse at the start of the year before natural mortality and growth occur. We believe this is a reasonable assumption given that the majority of natural mortality and growth of eastern oysters have been observed during spring and summer months (Shumway, 1996; Vølstad et al., 2008) The probability of growth and the natural mortality for smalls and markets each year were an estimated parameters. The natural mortality of spat was assumed to be known and constant at an instantaneous rate of 0.7 per year (50% survival) based on estimates of natural mortality from oysters in the Great Wicomico River, VA (Southworth et al., 2011). Similarly, the number of markets each year was calculated as the sum of the number of smalls that survive harvest and natural mortality and grow into markets and the number of markets already the previous year that survived natural mortality and harvest (Eq. 4), 𝑁𝑦+1,2 = (𝑁𝑦,2 − 𝐶𝑦,2)𝑒 −𝑀𝑦 + (𝑁𝑦,1 − 𝐶𝑦,1)𝐺𝑒 −𝑀𝑦. (4) Natural mortality for smalls and markets was estimated for each year as model parameters. The exploitation rate was calculated as the number of markets harvested divided by the estimated abundance of markets at the beginning of the year (Eq. 5), 𝑢𝑦 = 𝐶𝑦,2 𝑁𝑦,2 . (5) 14 The number of boxes each year was calculated by adding the number of boxes that did not disarticulate and were not broken by fishing activities to the number of adult oysters (small or market) that survived harvest but experienced natural mortality (Eq. 6), 𝐵𝑦+1,𝑠 = (𝑁𝑦,𝑠 − 𝐶𝑦,𝑠)(1 − 𝑒 −𝑀𝑦) + 𝐵𝑦,𝑠𝑒 −𝑏𝑦,𝑠(1 − 𝑢𝑦). (6) The model assumed that all smalls and markets became boxes after experiencing natural mortality. The model estimated the abundance of small and market boxes in the first year as parameters. Figure 3. Diagram of the fate of oysters in each stage of the population model. The relative area of oyster habitat was estimated by forcing an exponential decline through the estimated amount of habitat in 2010 (Eq. 7), ℎ𝑦 = 𝑒 −𝑑(𝑦−2010). (7) Oyster habitat in Chesapeake Bay has been modeled as an exponential decline previously (Wilberg et al. 2011), and an exponential decline model is appropriate given the degradation of oyster habitat that has been documented over time (Rothschild et al. 1994; Smith et al. 2005). The amount of restored habitat, was added to the estimated amount habitat each year (Eq. 8), 𝐻𝑦 = ℎ𝑦𝐻2011 + 𝐻𝑦,𝑎. (8) 15 2. Observation model Catchability was calculated for each of the five stages from the fall dredge survey and for the fishery-dependent indices of density. Catchability was calculated using a geometric mean difference between the observed index and estimated density on the log scale for each stage of the fall dredge survey (Eq. 9) and CPUE (Eq. 10), 𝑙𝑜𝑔𝑒(𝑞𝑋𝑠) = ∑ 𝑙𝑜𝑔𝑒(𝐼𝑋 𝑦,𝑠)−𝑙𝑜𝑔𝑒( 𝑋𝑦,𝑠 𝐻𝑦 )𝑦 𝑌 , (9) 𝑙𝑜𝑔𝑒(𝑞𝑔) = ∑ 𝑙𝑜𝑔𝑒(𝐼𝐶𝑏 𝑦)𝑦 −𝑙𝑜𝑔𝑒( 𝑁𝑦,2 𝐻𝑦 ) 𝑌 ., (10) which is the analytic solution for the maximum likelihood estimate. Catchability for the fishery dependent indices at the end of the fishing season, however, required adjusting abundance for the amount of exploitation (Eq. 11), 𝑙𝑜𝑔𝑒(𝑞𝑔) = ∑ 𝑙𝑜𝑔𝑒(𝐼𝐶𝑒 𝑦)𝑦 −𝑙𝑜𝑔𝑒( 𝑁𝑦,2(1−𝑢𝑦) 𝐻𝑦 ) 𝑌 . (11) The model predicted indices of density for the fall dredge survey and fishery CPUE. Predicted indices of density were calculated as the product of catchability and density at the beginning of the season for live oyster stages (Eq. 12), boxes (Eq. 13) and fishery CPUE (Eq. 14), 𝐼𝑁 𝑦,𝑠 = 𝑞𝑁,𝑠 𝑁𝑦,𝑠 𝐻𝑦 , (12) 𝐼𝐵 𝑦,𝑠 = 𝑞𝐵,𝑠 𝐵𝑦,𝑠 𝐻𝑦 , (13) 𝐼𝐶𝑏 𝑦,𝑔 = 𝑞𝑔 𝑁𝑦,2 𝐻𝑦 . (14) Predicted fishery CPUE at the end of the season required reducing the abundance by the amount of exploitation (Eq. 15), 16 𝐼𝐶𝑒 𝑦,𝑔 = 𝑞𝑔 𝑁 𝑦,2(1−𝑢) 𝐻𝑦 . (15) 3. Model fitting Model parameters were by minimizing the negative log likelihood function and penalties for some of the parameters. The objective function was the sum of the likelihood components and penalties (Eq. 16), −𝐿𝐿 = 𝐿𝑠𝑝 + 𝐿𝑠𝑚 + 𝐿𝑚𝑘 + 𝐿𝑠𝑚𝑏 + 𝐿𝑚𝑘𝑏 + 𝐿𝑎𝑏𝑑 + 𝐺𝑝 + 𝑏𝑝 + ?̌?𝑝 + ?̌?𝑝 + 𝑞𝑝 (16) We used a lognormal likelihood function (with additive constants ignored) for all indices in the model (Eq. 17), 𝐿𝑋 = 𝑙𝑜𝑔𝑒(𝜎𝑋) + ∑ 1 2 ( 𝑙𝑜𝑔𝑒(𝑋)−𝑙𝑜𝑔𝑒(?̂?) 𝜎𝑋 ) 2 𝑦 . (17) If data were not available for a year, that year was not included in the likelihood function. The models for Broad Creek, the Little Choptank River Closed, and Harris Creek Closed were fitted to estimates of abundance from patent tong monitoring in 2011 for Harris Creek and the Little Choptank, in 2013 for Broad Creek, and from the 2015 sanctuary monitoring for Harris Creek using a lognormal distribution (Eq. 18), 𝐿𝑎𝑏𝑑 = 1 2 ( 𝑙𝑜𝑔𝑒(𝑁𝑦𝑟,𝑎𝑏𝑠)−𝑙𝑜𝑔𝑒(𝑁𝑠𝑢𝑚−𝐶𝑠𝑢𝑚)𝑒 −2/3𝑀⃛ 𝜎𝑎𝑏𝑑 ) 2 . (18) The log-scale standard deviation was specified as 0.13 (the estimated standard error of abundance estimates from bootstrapping). Estimated abundance was calculated as the sum of live oysters minus oysters caught that year, and multiplied by the fraction of mortality that occurred between the end of the fishery and time of monitoring. Models assumed 2/3 natural mortality occurred between the end of the fishery and data collection. Accounting for the fraction of mortality was necessary because the majority of patent tong monitoring data were from late 17 August; by late August, roughly four of six months between the end of the previous fishery season and beginning of the next fishery season have elapsed. We incorporated penalties on some of the parameters to stabilize the estimates and to include outside information in the parameter estimation (Maunder, 2003). The probability of growth was penalized using a log-normal prior with a median of 0.45 based on growth of known age oysters in Maryland oyster sanctuaries (Paynter et al., 2010; Wilberg et al., 2011), 𝐺𝑝 = 1 2 ( 𝑙𝑜𝑔𝑒(?̅?)−𝑙𝑜𝑔𝑒(𝐺) 𝜎𝐺 ) 2 . (19) Lognormal penalties were also applied to the disarticulation rate of boxes (Eq. 20;(Christmas et al., 1997)), 𝑏𝑝 = 1 2 ( 𝑙𝑜𝑔𝑒(?̅?)−𝑙𝑜𝑔𝑒(𝑏) 𝜎𝑏 ) 2 , (20) deviations from mean mortality on the log scale (Eq. 21), ?̌?𝑝 = ∑ 1 2 ( 𝑙𝑜𝑔𝑒 ?̌?𝑦 𝜎?̌? ) 2 𝑦 , (21) and deviations from mean recruitment (Eq. 22), ?̌?𝑝 = ∑ 1 2 ( 𝑙𝑜𝑔𝑒 ?̌?𝑦 𝜎?̌? ) 2 𝑦 . (22) Catchability coefficients for each stage of the fall dredge survey were penalized with lognormal priors using the Broad Creek model estimates of catchability, 𝑞𝑝 = 1 2 ( 𝑙𝑜𝑔𝑒(𝑞)−𝑙𝑜𝑔𝑒(𝑞𝐵𝐶) 𝜎𝑞 ) 2 , (23) because the unconstrained models tended to estimate population size and rates of exploitation that were inconsistent with the substantial reduction in CPUE during the fishing season. The log-scale standard deviations (SD) for each penalty were assumed to be known. The log-scale SDs for the disarticulation rate of boxes were assumed to be 0.7. The log-scale SD for 18 growth was assumed to be 0.3 following Wilberg et al. (2011); the log-scale SD for the natural mortality deviations was assumed to be 0.5 (Shumway, 1996). The log-scale SD for deviations in recruitment was assumed to be 2.0 due to the highly variable nature of the recruitment process (North et al., 2008; Southworth et al., 2010). Table 1. Definition of model variables and parameters. Parameter Definition Indicator variables and subscripts y Year s Stage bar Oyster bar used as random effect in index standardization 0 Spat 1 Small 2 Market a Added with respect to habitat, spat p Subscript denoting prior g Gear G Growth subscript for prior abd Abundance subscript for absolute abundance likelihood component calculation abs Subscript for absolute abundance sum Subscript for sum of abundance and harvest quantities X Subscript denoting stage calculated in likelihood function Stage subscripts sp Subscript denoting spat oysters sm Subscript denoting small oysters mk Subscript denoting market oysters smb Subscript denoting small boxes mkb Subscript denoting market boxes Estimated parameters ?̅? Mean recruitment M Natural mortality of small and market oysters ?̅? Average rate of natural mortality for small and market oysters b Rate of disarticulation of small boxes and market boxes G Mean growth parameter d Rate of habitat decay 𝑑𝑠𝑚 Rate of disarticulation of small boxes 𝑑𝑚𝑘 Rate of disarticulation of market boxes 𝑁0 Initial abundance of small and market oysters 19 𝐵0 Initial abundance of small boxes and market boxes 𝐼𝑁 Estimated index of relative density of spat, small, and market oysters 𝐼𝐵 Estimated index of relative density of small boxes and market boxes 𝐼𝐶𝑏 Estimated index of CPUE at the beginning of the fishery 𝐼𝐶𝑒 Estimated index of CPUE at the end of the fishery Calculated quantities N Abundance of spat, small, and market oysters ?̅? Average abundance of small and market oysters B Abundance of small boxes and market boxes h Relative habitat H Absolute habitat u Rate of exploitation 𝑞𝑁 Catchability of spat, small, and market oysters 𝑞𝐵 Catchability of small boxes and market boxes 𝑞𝑔 Catchability by gear: hand tong, power dredge L Negative log-likelihood equation -LL Negative log-likelihood objective function Data C Harvest of market oysters, some small oysters S Rate of survival for planted spat 𝐼𝑁 Index of relative density for spat, small, and market oysters 𝐼𝐵 Index of relative density for small boxes and market boxes 𝐼𝐶𝑏 Index of catch per unit effort (CPUE) at beginning of fishery 𝐼𝐶𝑒 Index of catch per unit effort (CPUE) at end of fishery 𝑀𝑠𝑝 Natural mortality of spat 𝑁𝑎 Added spat on shell 𝑆𝑎 Survival of added spat on shell 𝐻𝑎 Added habitat ?̌? Annual recruitment deviations ?̌? Annual deviations in natural mortality for small and market oysters ?⃛? Fraction of natural mortality that occurred before patent tong survey 𝜎𝑋 Log-scale standard deviation for spat, small, market, small box, market box likelihood 𝜎?̌? Log-scale standard deviation of recruitment deviations 𝜎?̌? Log-scale standard deviation of natural mortality deviations 𝜎𝑏 Log-scale standard deviation of rate of box disarticulation 𝜎𝑞 Log-scale standard deviation of catchability 𝜎𝑎𝑏𝑑 Log-scale standard deviation of absolute abundance 𝜎𝐺 Log-scale standard deviation of growth probability 20 Sensitivity analyses We performed sensitivity analyses to determine how assumed input values affected model estimates. To assess the sensitivity of the model to the fraction of small oysters in the harvest, we ran the model with 2%, and 10% of harvest consisting of smalls. We also ran the model using different average sizes for smalls, 56 mm and 71 mm, that changed the number of smalls per bushel to 410 and 313, respectively. Additionally, we conducted two sensitivity runs with spat mortality at 0.3 and 0.7 per year. Lastly, we conducted several sensitivity analyses to determine the effects of the penalties on our estimates; these included doubling the SDs for growth, natural mortality deviations, and recruitment deviations. We ran sensitivity analyses using the Broad Creek model to avoid excessive model runs and because the Broad Creek model results were the most stable. We characterized the sensitivity of estimates of the probability of growth, the mean annual rate of natural mortality, the rate of habitat decline, mean recruitment, and the mean abundance of adults. Correlation tests Pearson correlation tests were run using R statistical software to determine if there were regional correlations in patterns of natural mortality and recruitment. Given that the stock assessments were conducted for regions within one large tributary, we were interested whether correlations in processes such as natural mortality and recruitment were observable among regions. 21 Results Model fits The model fit the pattern of the observed fall dredge survey indices well for all stages across all regions (Figures 4-8). Estimated relative density of spat followed the pattern of higher values at the beginning and ending periods of the time series, but were less variable than the observed values (Figure 4). Estimated relative density of both adult stages (live small and market sized oysters) showed a large decline in the early or mid-2000s and generally increased after 2005, which matched the trends in the data (Figures 5 & 6). The observed indices of small and market boxes were highest in the early 1990s and early 2000s (Figures 7 & 8) and were generally lower after 2004. Model estimates generally tracked the small and market box relative densities well, but the fit to small boxes was better than for market boxes. The estimates of both categories showed some serial over- and under-estimation in some regions. Figure 4. Estimated (line) and observed (points) log relative density of spat in number of oysters per bushel of cultch from the MD DNR fall dredge survey in all regions of the Choptank River Complex during 1989-2015. 22 Figure 5. Estimated (line) log relative density of small oysters in number of oysters per bushel of cultch fit to observed density (points) from the MD DNR fall dredge survey in all regions of the Choptank River Complex from 1989-2015. Figure 6. Estimated (line) log relative density of market oysters in number of oysters per bushel of cultch fit to observed density (points) from the MD DNR fall dredge survey in all regions of the Choptank River Complex during 1989-2015. 23 Figure 7. Estimated (line) log relative density of small boxes in number of boxes per bushel of cultch fit to observed density (points) from the MD DNR fall dredge survey in all regions of the Choptank River Complex during 1989- 2015. Figure 8. Estimated (line) log relative density of market boxes in number of boxes per bushel of cultch fit to observed density (points) from the MD DNR fall dredge survey in all regions of the Choptank River Complex during 1989-2015. 24 Estimated CPUE for hand tong and power dredge fisheries fit the observed values reasonably well, although observed values were not available for some years and regions. The trend in hand tong CPUE at the beginning (October) and end (March) of the fishing season (Figures 9 & 10) tracked with the trend in small and market relative density for each of the nine regions. Trends in the observed hand tong CPUE were similar for the beginning and end of the fishery, but the CPUE tended to be higher at the beginning than at the end. Power dredge CPUE was substantially greater at the beginning of the fishery in November (Figure 11) than at the end in March (Figure 12) in regions where this method of fishing was allowed. Figure 9. Estimated (line) and observed (points) of log catch per unit effort (CPUE) for the hand tong fishery at the beginning of the season (October) from the MD DNR buy ticket data for all regions of the Choptank River Complex during 1989-2015. 25 Figure 10. Estimated (line) and observed (points) of log catch per unit effort (CPUE) for the hand tong fishery at the end of the season (March) from the MD DNR buy ticket data for all regions of the Choptank River Complex during 1989-2015. Figure 11. Estimated (line) and observed (points) of log catch per unit effort (CPUE) for the power dredge fishery at the beginning of the season (November) from the MD DNR buy ticket data for all regions of the Choptank River Complex during 1989-2015 26 Figure 12. Estimated (line) and observed (points) of log catch per unit effort (CPUE) for the power dredge fishery at the end of the season (March) from the MD DNR buy ticket data for all regions of the Choptank River Complex during 1989-2015. Mortality The estimated natural mortality rate showed a consistent pattern among regions with mortality rates highest (50% - 90% per year) during the early 1990s and early 2000s (Figure 13). Natural mortality was lower during 2004 – 2015 and varied by region ranging from 6% to 27% per year, except for the Harris Creek Closed model; the Harris Creek Closed model estimated an average of 41% per year during 2004 – 2015, which was heavily influenced by the high estimate in the terminal year. Overall, the Harris Creek Closed had the highest average annual rate of natural mortality at 56%, and the Little Choptank Open model estimated the lowest average rate of natural mortality at 16%. 27 Figure 13. Estimated natural mortality (% per year) for small and market adult oysters and 95% confidence intervals (shaded) in each region of the Choptank River Complex during 1989 - 2015. Recruitment Estimated recruitment decreased from an average of 93 million spat in 1989 to 26 million spat in 2015, a 71% reduction (Figure 14). Broad Creek and the Harris Creek Closed regions had the highest average recruitment of 123 million and 90 million spat, respectively, across the time series. The Tred Avon and Upper Choptank models had the lowest average recruitment at 12 million and 23 million, respectively. In most regions, the highest estimated recruitment occurred in 1990 or in 1997, except for the Upper Choptank. Estimated recruitment was generally low during 2000-2010, after which a small increase occurred. For Broad Creek, the Lower Choptank, the Middle Choptank, the Upper Choptank, Little Choptank Closed, Harris Creek Open and Harris Creek Closed regions, recruitment increased in 2012. Estimated recruitment also 28 increased in the last year of model for the Lower Choptank, the Middle Choptank and Little Choptank Open regions. Figure 14. Estimated recruitment with 95% confidence intervals (shaded) in each region of the Choptank River Complex during 1989 – 2015. Abundance Estimated adult abundance, consisting of small and market oysters, was highest in the late 1980s or early 1990s for all of the regions except Broad Creek; the estimated peak abundance occurred in the late 1990s in Broad Creek (Figure 15). Total estimated adult abundance declined in the CRC by 50% during 1989-2015. Broad Creek and the Little Choptank Closed regions had the highest average abundance at 140 million, and 65 million adults, respectively, during 1989-2015. The Tred Avon and Harris Creek Open regions both had the lowest average adult abundance at 29 million and 44 million, respectively. Abundance increased in all regions in the late 1990s and then declined substantially thereafter with abundance in most regions being below their average at the end of the time series. The main exception was Broad 29 Creek, which experienced an increase in abundance starting in 2005. The other exception was the Harris Creek Closed region that increased by 110 million adult oysters after 2010. The increases in estimated adult abundance in these regions were preceded by high recruitment in the previous year. Figure 15. Estimated adult (small and market) eastern oyster abundance with 95% confidence intervals (shaded) in each region of the Choptank River Complex during 1989 – 2015. Habitat Patterns in habitat change differed among regions (Figure 16), with some regions suffering substantial amounts of habitat loss and others experiencing low amounts. On average, the amount of available oyster habitat in the CRC declined by an estimated 67% during 1989- 2015. The amount of estimated habitat declined the least in the Upper Choptank, Broad Creek, and the Middle Choptank. The steepest rates of habitat loss were in the Tred Avon (9% per year), 30 Lower Choptank (7% per year), and Harris Creek Closed (6 % per year). The Upper Choptank had the lowest rate of habitat loss (<0.01%). Figure 16. Estimated change in the area of hard bottom oyster habitat and 95% confidence intervals for each region of the Choptank River Complex during 1989 - 2015. Exploitation Estimated exploitation rates were highest during the same periods when adult abundance was highest for most regions (Figure 17). During 1989 – 2015 Harris Creek Closed experienced the highest average exploitation rate (14.8% per year) prior to becoming an oyster sanctuary, and the Tred Avon region experienced the lowest average exploitation rate (4.2% per year). The dominant pattern in exploitation showed moderate to high rates between 40% and 50% per year in the early 1990s. Exploitation increased from low rates to between 10% and 50% per year during 1995 – 2000. Exploitation rates below 10% were estimated after this period followed by another sharp increase after 2010 in most regions. Notable exceptions to this pattern included Broad Creek, where estimated exploitation rate was higher during 2005-2010 and 2012 – 2015 31 (20-45% per year), Harris Creek Closed, where the exploitation was highest during 1995 – 2000 and 2005 – 2010 (50-60% per year), and the Lower Choptank and Harris Creek Open, where the exploitation rate increased substantially after 2010. The estimated exploitation rate remained low in the Middle and Upper Choptank, Tred Avon, and Little Choptank Closed. The average annual rate of exploitation among all regions of the CRC during 1989 – 2015 was 9.1%. For most years, exploitation was a smaller proportion of the total mortality (<0.5) than natural mortality. However, during years with the highest abundance of adult oysters, the proportion of mortality due to exploitation was greater than natural mortality (>0.5); this was true for all regions except for the Tred Avon (Figure 18). Figure 17. Estimated exploitation rate (% per year) of market-sized eastern oysters with 95% confidence intervals (shaded) for all regions of the Choptank River Complex during 1989 - 2015. 32 Figure 18. Proportion of mortality due to exploitation (e.g., 0 indicates no mortality from fishing, 0.5 indicates 50% of mortality was due to fishing). Parameter estimates Estimates of parameter values were similar among regions for most parameters (Table 2). The probability of growth from the small to market stage was consistently between 0.30 and 0.50 per year for all regions except for the Lower Choptank, which had the lowest value (0.26 per year). The catchability for live small and market stages was generally 2-3 times higher than the catchability of boxes or spat. Estimates of catchability for the Lower Choptank model were higher than the other regions. Instantaneous rates of disarticulation were consistently between 0.8 and 1.5 per year and were greater for small boxes than market boxes. The highest rates of disarticulation were estimated in the Lower Choptank model at 1.6 per year for small boxes, and 1.3 per year for market boxes. 33 Table 2. Estimated values of parameters that were assumed constant in the model. G represents the probability of an oyster growing from the small to market stage (per year), h the rate of habitat loss (% per year), d the instantaneous rate of box disarticulation (per year), and q the catchability coefficient. The subscripts for d and q indicate the stage: spat (sp), small (sm), market (mk), small box (smb), and market box (mkb). Regions are identified using the following abbreviations: Broad Creek (BC), Lower Choptank River (LC), Middle Choptank River (MC), Upper Choptank River (UC), Tred Avon River (TA), Little Choptank Open (LCO), Little Choptank Closed (LCC), Harris Creek Open (HCO), and Harris Creek Closed (HCC). The subscripts for d and q indicate the stage: spat (sp), small (sm), market (mk), small box (smb), and market box (mkb). G h (%) 𝑑𝑠𝑚 𝑑𝑚𝑘 𝑞𝑠𝑝 𝑞𝑠𝑚 𝑞𝑚𝑘 𝑞𝑠𝑚𝑏 𝑞𝑚𝑘𝑏 BC 0.30 1.25 1.03 0.85 5.30 15.05 10.94 3.84 4.48 LC 0.26 7.14 1.59 1.25 13.92 26.01 42.49 7.53 21.02 MC 0.50 2.01 0.87 0.83 3.31 13.52 9.79 4.26 7.73 UC 0.28 0.00 1.24 1.31 2.40 4.78 4.06 1.47 2.88 TA 0.48 9.24 0.98 0.85 3.14 14.61 10.37 4.59 5.35 LCO 0.48 2.53 0.85 0.74 3.28 13.30 9.40 4.61 8.28 LCC 0.38 2.41 1.35 1.37 3.08 13.19 11.68 3.83 6.88 HCO 0.35 5.26 1.27 1.00 4.13 16.69 13.76 3.45 4.38 HCC 0.31 5.86 0.98 1.35 3.03 17.57 40.52 3.66 12.31 Sensitivity analyses The Broad Creek model was not sensitive to alternative assumptions about the fraction of small oysters in the harvest, the average size of small oysters in the harvest, the natural mortality rate for spat, or the standard deviations for natural mortality, recruitment variability, or growth (Table 3). The parameter estimates from the sensitivity scenarios were almost always within 30% of the estimates for the base model. Model estimates were moderately sensitive in the scenario with higher variability of natural mortality; the average natural mortality rate, average recruitment, and average abundance decreased by approximately 30%. Generally, mean recruitment and adult abundance were the most sensitive estimates. The estimate of habitat loss for Broad Creek was relatively low (1.25% per year), so even large percentage-wise changes in the estimated rate of habitat loss were on the order of +/- 2% per year differences. 34 Table 3. Sensitivity analyses results for the Broad Creek estimation model. Values indicate the relative differences (%) from the parameter values obtained from the base model. Parameter values are from left to right: growth (G), average annual rate of natural mortality (M), average annual recruitment (r), and average annual adult abundance (N), and rate of habitat loss (h). Please, refer to the methods section for scenario descriptions. Scenario G M r N h 2% smalls 0.41 0.22 2.55 2.54 2.82 10% smalls -0.81 -0.66 -4.19 -4.18 -4.97 56 mm smalls -0.13 -0.10 0.26 0.21 0.30 71 mm smalls 0.05 -0.11 -0.20 -0.16 -0.21 𝑀𝑠𝑝 = 0.5 -0.39 -1.56 -19.28 -0.49 -16.75 𝑀𝑠𝑝 = 0.9 0.28 1.38 24.16 0.72 15.86 𝜎?̌? (2x) 14.00 -42.28 -47.92 -36.08 -99.64 𝜎?̌? (2x) -3.50 5.51 25.96 21.03 183.55 𝜎𝐺 (2x) -24.49 7.19 32.41 26.82 11.18 Correlation tests Correlation tests indicated similar patterns of natural mortality and recruitment among regions. All correlations for natural mortality were significant and positive (Table 4). Most correlations for recruitment were also significant and positive, with the strongest correlations in recruitment between geographically adjacent regions (Table 5; e.g. Broad Creek and the Lower Choptank, Lower Choptank and Middle Choptank, Little Choptank Open and Closed, Harris Creek Open and Closed). 35 Table 4. Pearson correlation coefficients (above the diagonal) and p-values (shaded region below the diagonal) for the natural mortality rate among regions. Regions are identified using abbreviations: Broad Creek (BC), Lower Choptank River (LC), Middle Choptank River (MC), Upper Choptank River (UC), Tred Avon River (TA), Little Choptank Open (LCO), Little Choptank Closed (LCC), Harris Creek Open (HCO), and Harris Creek Closed (HCC). BC LC MC UC TA LCO LCC HCO HCC BC 1.00 0.60 0.67 0.47 0.85 0.60 0.72 0.70 0.73 LC 0.01 1.00 0.89 0.70 0.76 0.80 0.92 0.95 0.55 MC <0.01 <0.01 1.00 0.65 0.87 0.96 0.91 0.79 0.59 UC 0.04 <0.01 <0.01 1.00 0.50 0.56 0.65 0.66 0.44 TA <0.01 <0.01 <0.01 0.04 1.00 0.87 0.85 0.74 0.53 LCO 0.01 <0.01 <0.01 0.02 <0.01 1.00 0.86 0.69 0.41 LCC <0.01 <0.01 <0.01 <0.01 <0.01 <0.01 1.00 0.91 0.60 HCO <0.01 <0.01 <0.01 <0.01 <0.01 <0.01 <0.01 1.00 0.61 HCC <0.01 0.02 0.01 0.05 0.03 0.05 0.01 0.01 1.00 36 Table 5. Pearson correlation coefficients (above the diagonal) and p-values (shaded region below the diagonal) for recruitment among regions. Regions are identified using abbreviations: Broad Creek (BC), Lower Choptank River (LC), Middle Choptank River (MC), Upper Choptank River (UC), Tred Avon River (TA), Little Choptank Open (LCO), Little Choptank Closed (LCC), Harris Creek Open (HCO), and Harris Creek Closed (HCC). BC LC MC UC TA LCO LCC HCO HCC BC 1.00 0.57 0.42 -0.19 0.64 0.34 0.65 0.77 0.36 LC 0.03 1.00 0.95 -0.02 0.90 0.89 0.81 0.89 0.47 MC 0.36 <0.01 1.00 -0.08 0.83 0.98 0.79 0.82 0.35 UC 1.00 1.00 1.00 1.00 -0.03 -0.11 -0.11 -0.11 0.07 TA <0.01 <0.01 <0.01 1.00 1.00 0.75 0.85 0.92 0.71 LCO 0.84 <0.01 <0.01 1.00 <0.01 1.00 0.77 0.74 0.32 LCC <0.01 <0.01 <0.01 1.00 <0.01 <0.01 1.00 0.86 0.56 HCO <0.01 <0.01 <0.01 1.00 <0.01 <0.01 <0.01 1.00 0.59 HCC 0.74 0.18 0.78 1.00 <0.01 0.96 0.04 0.02 1.00 Discussion Abundance of eastern oyster in the CRC has increased since 2004 because of relatively low natural mortality rates and increased recruitment in recent years. Natural mortality rates of eastern oysters have been stable below 30% per year in almost all regions of the CRC since 2004, the longest period of low mortality in 27 years (Tarnowski, 2016). This was potentially due to a confluence of advantageous environmental conditions, and an increased resistance to Dermo (Kimmel & Newell, 2007; North et al., 2008; Yu et al., 2011; Zhang et al., 2014). Recruitment events may have been influenced by salinity (Kimmel & Newell, 2007) and freshwater flow (Tarnowski, 2016), and were characterized large recruitment events in the late 1990s, 2010, and 2012. The high recruitment events caused increased abundance during 2010 – 2015 in regions of the CRC open to harvest. Exploitation rates were generally above sustainable 37 rates (5-10%; Wilberg et al., 2013) during 2013 – 2015 in the regions where fishing is allowed, which suggests that the oyster recovery may not be sustained. The pattern of natural mortality rates over time matched closely with average mortality rate estimates across oyster bars in Chesapeake Bay, Maryland, from Jordan and Coakley (2004), Vølstad et al. (2008), and Wilberg et al. (2011). Previous studies estimated a lower natural mortality rate, 50%, per year during 2000 – 2002 (Jordan & Coakley, 2004; Vølstad et al., 2008; Wilberg et al., 2011), whereas we estimated rates between 50 and 90% per year during that same period. Higher natural mortality rates in the CRC than the Maryland Chesapeake Bay, on average, are consistent with intensity and prevalence of Dermo from disease bar sentinel site observations (Tarnowski, 2016). Additionally, Vølstad et al. (2008) and Wilberg et al. (2011) estimated consistent, lower rates of natural mortality around 25% after 2004; our results showed that low natural mortality rates persisted with an average rate of annual mortality between 6 – 27% during 2004 – 2015 for the CRC as a whole, excluding Harris Creek Closed region. The Harris Creek Closed model’s 41% average natural mortality was affected by high estimates in the two most recent years. These high mortality rates may have been caused by fitting the model to estimated abundance in 2014 from patent tong surveys to monitor the progress of restoration efforts. In addition, our results may provide a more accurate estimate of natural mortality than estimates from other approaches. The way natural mortality was modeled within our study differed from that of Vølstad et al. (2008); our models estimated different catchabilities between live and box stages of oysters similar to Powell et al., (2007). Additionally, our method did not assume that all boxes disarticulate within one year’s time (Vølstad et al., 2008). These changes made the estimation of boxes in each region’s population more consistent with the observed fall dredge indices, and did not count all boxes as recent deaths. 38 A combination of advantageous environmental conditions and the development of resistance to Dermo may be responsible for low natural mortality since 2004. Salinity and freshwater flow have been near their average during 2003 – 2015, Chesapeake Bay-wide, but Dermo has remained prevalent throughout the CRC, albeit at relatively low intensity (Tarnowski, 2016). The strain of Dermo in Chesapeake Bay may be particularly lethal (Brown, et al., 2005). However, exploration of the eastern oyster genetics has suggested that eastern oysters possess multiple innate mechanisms to combat the physiological effects of Dermo (Yu et al., 2011; Zhang et al., 2014). The cvSI-1 gene in eastern oyster has a highly polymorphic response to Dermo that inhibits Dermo’s ability to proliferate and characterizes the Dermo-resistant strain of eastern oyster (Yu et al., 2011). Eastern oysters also possess a high degree of genetic diversity with respect to combating apoptosis (cell death), another symptom of Dermo (Zhang et al., 2014). Years with high recruitment were fairly consistent among regions of the CRC. Large recruitment events were evident in the late 1990s, 2010, and 2012 in most regions. These years also had high recruitment in other regions of Maryland (Tarnowski, 2016). Inter-annual variability in recruitment is thought to be heavily influenced by environmental factors, primarily by salinity regime (Kimmel & Newell, 2007); salinity influences oyster growth (Shumway, 1996) and larval swimming ability, which can also improve the success of settlement (North et al., 2008). These patterns likely affected all regions of the CRC and may explain the strong correlation in recruitment among regions. The CRC is similar to other Chesapeake Bay tributaries like the Great Wicomico River, VA in that recruitment events are episodic, and drive population demographics (Southworth et al., 2010). Like Delaware Bay, New Jersey, strong recruitment events in the CRC tend to be associated with low freshwater flow events (North et 39 al., 2008; Narvaez et al., 2012; Tarnowski, 2016), and high salinity during summer months conducive to successful recruitment (Kimmel & Newell, 2007; 2012). Our estimates of abundance in the CRC were consistent with previous Maryland-wide estimates of abundance and the fraction of harvest that came from the CRC (Tarnowski 2016). Wilberg et al. (2011) estimated the abundance of adult (small and market stages) eastern oysters in the Maryland portion of Chesapeake Bay to be 851 million in 2009; our models estimated abundance after 2009 to be between 30% - 40% of that value, which was close to the average fraction of harvest from the CRC (28%; Tarnowski 2016). Most of from the abundance in the CRC was in Broad Creek with an average of 160 million adults after 2010, which supports 19% of the eastern oyster harvest from the Maryland portion of the Chesapeake Bay (Tarnowski, 2016). Both models use harvest information to inform population size, but our models also included more data than was used in Wilberg et al. (2011). Patterns in exploitation have changed most over 2010 – 2015 and were likely due to a combination of factors including implementation of large sanctuaries (i.e., marine protected areas), increased power dredging effort, and strong recruitment events in 2010 and 2012. The average exploitation rate was 9.1% among regions for the time series, but patterns over time differed among regions. The rate of exploitation in Broad Creek increased from 2.9% to 31.6%, the Lower Choptank increased from 2.7% to 71%, and Harris Creek Open increased from 3.2% to 25.2% during 2010 - 2015. These increased rates were associated with substantial increases in effort; during 2010 - 2015, power dredging effort in the Lower Choptank and Harris Creek Open increased ten-fold, and four-fold, respectively. Broad Creek appears to have experienced a greater rate of exploitation than other CRC regions during 2010 – 2015. The rates of exploitation in the CRC during 2010-2015 were generally 20% per year or higher during years of high 40 abundance, which was about twice the upper limit of the exploitation rate that would maximum sustainable yield (10%; Wilberg et al. 2013), and much higher than is allowed in Delaware Bay, 2% per year (Bushek & Ashton-Alcox, 2013). Thus, exploitation rates in recent years may be above sustainable levels. The rate at which habitat declined was region-specific, and showed the most uncertainty of estimated quantities. Estimates of the annual rate of habitat decline were highest in Tred Avon (9% per year) and Lower Choptank (7% per year) regions. The rates from these regions are consistent with the amount of habitat degradation reported by Smith et al. (2005) in the Lower Choptank, Middle Choptank, and Tred Avon during 1999 - 2001. Estimates from regions that were predominantly fished with hand tongs, like the Middle Choptank and Upper Choptank, and Broad Creek, had lower estimated habitat loss (1 – 3% per year) over the 27-year time series. Although Harris Creek Closed had an estimated rate of 5.8% habitat loss per year, this trend abated during 2012 – 2014 due to extensive habitat restoration efforts made to the region (Westby et al., 2016; NOAA, 2016; Tarnowski, 2016). Wilberg et al. (2011) estimated a rate of 4% habitat loss per year for the Maryland portion of Chesapeake Bay during 1980-2008, which was consistent with the average rate in the CRC, 3.97% per year. The high contrast in estimated rate of habitat decline among regions may highlight a lack of ability by our models to capture the complexity of processes affecting the decline. Our models incorporated more data sources than previous stock assessments on eastern oysters in Chesapeake Bay (Wilberg et al. 2011) and were conducted at a finer spatial scale. We included fishery CPUE, estimates of habitat data, amount of habitat added by restoration activities, estimates of abundance from patent tong monitoring, and the number of eastern oysters planted. For stock assessments at a finer spatial scale, these data likely provided models 41 with vital information to accurately estimate quantities of interest. Furthermore, we used an updated estimate of the number of oysters per bushel and incorporated harvest of smalls into the model. Fitting the models to fishery CPUE provided important information about how much eastern oyster abundance was depleted during the fishing season. Our models used two rates of underreporting as well, which are likely to be more accurate than the assumed rate of 50% from Wilberg et al. (2011). These rates did not affect trends in population dynamics, but helped provide estimates on a more realistic scale given conversations with the OysterFutures work group. We also believe that including the harvest of smalls and updating the number of oysters per bushel to scale harvest resulted in more realistic estimates of exploitation and abundance. Our study represents the first eastern oyster stock assessments for Maryland’s three largest oyster sanctuaries. Stock assessments that include no-take areas without consideration of spatial structure often result in biased estimates because the closed areas mean the whole population is not subject to the same fishing mortality (Punt & Methot, 2004; Field et al., 2006; Pincin & Wilberg, 2012). To address this problem, we developed separate models for regions containing sanctuaries. Our models for the Little Choptank and Harris Creek tributaries account for the spatial differences in harvest activity by assuming zero harvest after 2010 for the closed portions. Spatial disaggregation at a fine spatial scale can be difficult when the data are insufficient (Fulton et al., 2015). For example, the Middle Choptank, Upper Choptank, and Tred Avon regions also contain oyster sanctuaries, but disaggregated models could not converge on a solution for multiple parameters. Harvest has been low in these regions for most of the 27 years modeled, and the models used in this study rely on robust fishing data to characterize the population dynamics of eastern oysters (Wilberg et al., 2011). 42 Habitat, recruitment, and abundance have declined by 50 – 70% in the CRC during the last 27 years, despite increases in abundance in recent years. Rates of exploitation increased beyond a sustainable rates during years of high abundance, particularly in regions where power dredging is the dominant gear type for fishing. The rates of habitat decline generally agree with the Beck et al. (2011) estimate of the loss of oyster reefs, worldwide (2011). Low natural mortality in recent years may indicate increased resistance to Dermo, which could lead to future increases in abundance. The fishery in the CRC appears to depend on high recruitment events, such that increases in abundance seem to lead to an increase in the exploitation rate. This type of fishery response may prevent the population from increasing substantially in regions that are open to harvest. 43 Appendices Appendix 1. Table of spat on shell planting by the Oyster Recovery Partnership (ORP) by region. Regions are identified using abbreviations: Broad Creek (BC), Lower Choptank River (LC), Middle Choptank River (MC), Upper Choptank River (UC), Tred Avon River (TA), Little Choptank Open (LCO), Little Choptank Closed (LCC), Harris Creek Open (HCO), and Harris Creek Closed (HCC). Region Year of planting Spat on shell planted (millions) LC 1997 1.02 2000 18.50 LC Total Planting 19.52 MC 1997 1.02 2000 0.84 2001 1.94 2003 17.50 2008 30.20 MC Total Planting 51.50 UC 1997 2.50 1998 4.77 2001 18.92 2002 14.93 2003 79.90 2004 17.10 2005 44.52 2006 91.20 2007 57.00 2008 129.78 UC Total Planting 460.62 TA 2009 14.07 2015 10.18 TA Total Planting 24.25 LCC 2002 3.25 44 2003 2.30 2004 0.94 2011 45.66 2014 81.34 LCC Total Planting 127.00 HCC 2011 55.63 2012 450.00 2013 723.10 2014 444.30 2015 389.54 HCC Total Planting 2062.57 Total Planting by ORP in CRC 2725.94 45 References Allen, S., O’Neill, C., Sowers, A., Weissberger, E., & Westby, S. R. (2013). Harris Creek oyster restoration tributary plan: A blueprint to restore the oyster population in Harris Creek, a tributary of the Choptank River on Maryland’s Eastern Shore. Maryland Interagency Oyster Restoration Workgroup. Beck, M. W., Brumbaugh, R. D., Airoldi, L., Carranza, A., Coen, L. D., Crawford, C., … Guo, X. (2011). Oyster reefs at risk and recommendations for conservation, restoration, and management. Bioscience, 61(2), 107–116. Bocking, S. (2011). The oyster question: Scientists, watermen, and the Maryland Chesapeake Bay since 1880 by Keiner. Isis, 102(1), 147–148. Bolker, B. M., Brooks, M. E., Clark, C. J., Geange, S. W., Poulsen, J. R., Stevens, M. H. H., & White, J. S. (2009). Generalized linear mixed models: a practical guide for ecology and evolution. Ecology and Evolution, 24(3), 127–135. Brown, B. L., Butt, A. J., Shelton, S. W., Meritt, D., & Paynter, K. T. (2005). Resistance of Dermo in eastern oysters, Crassostrea virginica (Gmelin), of North Carolina but not Chesapeake Bay heritage. Aquaculture Research, 36, 1391–1399. Burreson, E. M., & Ragon Calvo, L. M. (1996). Epizootiology of Perkinsus marinus disease of oysters in Chesapeake Bay, with emphasis on data since 1985. Journal of Shellfish Research, 15(1), 17–34. Burreson, E. M., Stokes, N. A., & Friedman, C. S. (2000). Increased virulence in the introduced pathogen: Haplosporidium nelsoni (MSX) in the eastern oyster Crassostrea virginica. Journal of Aquatic Animal Health, 12(1), 1–8. 46 Bushek, D., & Ashton-Alcox, K. A. (2013). Stock assessment workshop New Jersey Delaware Bay oyster beds (15th SAW). Rutgers New Jersey Agricultural Experiment Station: Haskin Shellfish Research Laboratory. Chai, A., Homer, M., Tsai, C., & Goulletquer, P. (1992). Evaluation of oyster sampling efficiency of patent tongs and an oyster dredge. North American Journal of Fisheries Management, 12(4), 825–832. Christmas, J. F., McGinty, M. R., Randle, D. A., Smith, G. F., & Jordan, S. A. (1997). Oyster shell disarticulation in three Chesapeake Bay tributaries. Journal of Shellfish Research, 16(1), 115–123. Coen, L. D., Brumbaugh, R. D., Bushek, D., Grizzle, R., Luckenbach, M. W., Posey, M. H., … Tolley, S. G. (2007). Ecosystem services related to oyster restoration. Marine Ecology Progress Series, 341, 303–307. Coen, L. D., Luckenbach, M. W., & Breitburg, D. L. (1999). The role of oyster reefs as essential fish habitat: a review of current knowledge and some new perspectives. American Fisheries Society Symposium, 34, 303–307. Cook, T., Klinck, J., & Miller, J. (1998). The relationship between increasing sea-surface temperature and the northward spread of Perkinsus marinus (Dermo) disease epizootics in oysters. Estuarine, Coastal and Shelf Science, 46, 587–597. Field, J. C., Punt, A. E., Methot, R. D., & Thompson, C. J. (2006). Does MPA mean “major problem for assessments”? Considering the consequences of place-based management systems. Fish and Fisheries, 7, 284–302. Fournier, D., Skaug, H. J., Ancheta, J., Ianelli, J., Magnusson, A., Maunder, M. N., … Sibert, J. (2012). AD Model Builder: using automatic differentiation for statistical inference of 47 highly parameterized complex nonlinear models. Optimization Methods and Software, 27(2), 233–249. Fulton, E. A., Bax, N. J., Bustamante, R. H., Dambacher, J. M., Dichmont, C., Dunstan, P. K., … Smith, D. C. (2015). Modelling marine protected areas: insights and hurdles. Philosophical Transactions Royal Society B, 370. Retrieved from http://dx.doi.org/10.1098/rstb.2014.0278 Ihde, T. F., Wilberg, M. W., Secor, D. H., & Miller, T. J. (2011). Harnessing the knowledge of stakeholders to enhance U.S. marine recreational fisheries with to the Atlantic King Mackerel fishery. American Fisheries Society Symposium, 73, 75–93. Jordan, S. J., & Coakley, J. M. (2004). Long-term projections of eastern oyster populations under various management scenarios. Journal of Shellfish Research, 23(1), 63–72. Jordan, S. J., Greenhawk, K. N., McCollough, C. B., Vanisko, J., & Homer, M. (2002). Oyster biomass abundance, and harvest in northern Chesapeake Bay: trends and forecasts. Journal of Shellfish Research, 21, 733–741. Kennedy, V. S., & Breisch, L. L. (1983). Sixteen decades of political management of the oyster fishery in Maryland’s Chesapeake Bay. Journal of Environmental Management, 164, 153–171. Kimmel, D. G., & Newell, R. I. E. (2007). The influence of climate variation on eastern oyster (Crassostrea virginica) juvenile abundance in Chesapeake Bay. Limnological Oceanography, 52(3), 959–965. Lenihan, H. S., Micheli, F., Shelton, S. W., & Peterson, C. H. (1999). The influence of multiple stressors on susceptibility to parasites: An experimental determination with oysters. Limnological Oceanography, 44(3, part 2), 910–924. 48 Lenihan, H. S., & Peterson, C. H. (1998). How habitat degradation through fishery disturbance enhances impacts on hypoxia on oyster reefs. Ecological Applications, 8(1), 128–140. Linden, A., & Mantyniemi, S. (2011). Using the negative binomial distribution to model overdispersion in ecological count data. Ecology, 92(7), 1414–1421. Mann, R., & Powell, E. N. (2007). Why oyster restoration goals in the Chesapeake Bay are not and probably cannot be achieved. Journal of Shellfish Research, 26(4), 905–917. Maunder, M. N. (2003). Paradigm shifts in fisheries stock assessment: From integrated analysis to Bayesian analysis and back again. Natural Resource Modeling, 16(4), 465–475. Maunder, M. N., & Punt, A. E. (2004). Standardizing catch and effort data: a review of recent approaches. Fisheries Research, 70, 141–159. Miller, T. J., Blair, J. A., Ihde, T. F., Jones, R. M., Secor, D. H., & Wilberg, M. W. (2010). FishSmart: An innovative role for science in stakeholder-centered approaches to fisheries management. Fisheries, 35(9), 424–433. Narvaez, D. A., Klinck, J. M., Powell, E. N., Hofmann, E. E., Wilkin, J., & Haidvogel, D. B. (2012). Circulation and behavior controls on dispersal of eastern oyster (Crassostrea virginica) larvae in Delaware Bay. Journal of Marine Research, 70(2–3), 411–440. Newell, R. I. E. (1988). Ecological changes in Chesapeake Bay: Are they the result of overharvesting the American oyster, Crassostrea virginica? In Understanding the Estuary: Advances in Chesapeake Bay Research (pp. 536–546). Baltimore, MD. North, E. W., Schlag, Z., Hood, R. R., Li, M., Zhong, L., Gross, T., & Kennedy, V. S. (2008). Vertical swimming behavior influences the dispersal of simulated oyster larvae in a coupled particle-tracking and hydrodynamic model of Chesapeake Bay. Marine Ecology Progress Series, 359, 99–115. 49 Paynter, K. T., Michaelis, A., & Handschy, A. (2014). Paynter Lab Annual Monitoring and Research Summary 2013: Submitted to the Oyster Recovery Partnership. College Park, MD: University of Maryland Center for Environmental Science. Paynter, K. T., Politano, V., Lane, H. A., Allen, S. M., & Meritt, D. (2010). Growth rates and prevalence of Perkinsus marinus prevalence in restored oyster populations in Maryland. Journal of Shellfish Research, 29(2), 309–317. Pincin, J. S., & Wilberg, M. W. (2012). Surplus production model accuracy in populations affected by no-take marine protected area. Marine and Coastal Fisheries, 4(1), 511–525. Powell, E. N., Ashton-Alcox, K. A., & Kraeuter, J. N. (2007). Revaluation of oyster dredge efficiency in survey mode: Application in stock assessment. North American Journal of Fisheries Management, 27, 492–511. Punt, A. E., & Methot, R. D. (2004). Effects of MPAs on the assessment of marine fisheries. In Aquatic Protected Areas as Fisheries Management Tools (pp. 133–154). Quebec, Canada. Rothschild, B. J., Ault, J. S., Goulletquer, P., & Heral, M. (1994). Decline of the Chesapeake Bay oyster population: a century of habitat destruction and overfishing. Marine Ecology Progress Series, 111, 29–39. Shumway, S. E. (1996). Natural Environmental Factors. In The Eastern Oyster Crassostrea virginica (pp. 467–503). College Park, MD: Maryland Sea Grant. Smith, G. F., Bruce, D. G., Roach, E. B., Hansen, A., Newell, R. I. E., & McManus, A. E. (2005). Assessment of recent habitat conditions of eastern oyster Crassostrea virginica bars in mesohaline Chesapeake Bay. North American Journal of Fisheries Management, 25(4), 1569–1590. 50 Southworth, M., Harding, J. M., Wesson, J. A., & Mann, R. (2010). Oyster (Crassostrea virginica, Gmelin 1791) population dunamics on public reefs in the Great Wicomico River, Virginia USA. Journal of Shellfish Research, 29(2), 271–290. Tarnowski, M. T. (2016). Maryland oyster population status report. Annapolis, MD: Maryland Department of Natural Resources. Volstadt, J. H., Dew, J., & Tarnowski, M. T. (2008). Estimation of annual mortality rates for eastern oysters (Crassostrea virginica) in Chesapeake Bay based on box counts and application of those rates to project population growth of C. virginica and C. ariakensis. Journal of Shellfish Research, 27(3), 131–144. Westby, S. R., Bruce, D. B., Slacum, W., French, E., Sowers, A., & Weissberger, E. (2016). 2015 oyster restoration implementation update. Silver Spring, MD: NOAA. Wilberg, M. W., Livings, M. E., Barkman, J. S., Morris, B. T., & Robinson, J. M. (2011). Overfishing, disease, habitat loss, and potential extirpation of oysters in Upper Chesapeake Bay. Marine Ecology Progress Series, 436, 131–144. Wilberg, M. W., Wiedenmann, J. R., & Robinson, J. M. (2013). Sustainable exploitation and management of autogenic ecosystem engineers: application to oysters in Chesapeake Bay. Ecological Applications, 23(4), 766–776. Yu, H., He, Y., Wang, X., Zhang, Q., Bao, Z., & Guo, X. (2011). Polymorphism in a serine protease inhibitor gene and its association with disease resistance in the eastern oyster (Crassostrea virginica Gmelin). Fish and Shellfish Immunology, 1–6. Zhang, L., Li, L., Zhang, G., & Guo, X. (2014). Transcriptome analysis reveals a rich gene set related to innate immunity in the eastern oyster (Crassostrea virginica). Marine Biotechnology, 16, 17–33. 51