Analysis of fish population size distributions confirms cessation of fishing in marine protected areas

The number of protected areas that restrict or prohibit harvest of wild populations is growing. In general, protected areas are expected to increase the abundance of previously‐harvested species. Whether a protected area achieves this expectation is typically evaluated by assessing trends in abundance after implementation. However, the underlying assumption that harvest has actually ceased is rarely tested directly. Determining whether illegal harvest (poaching) has continued in a protected area is important to planning enforcement and adaptive management. Here, we estimated harvest rates for four kelp forest fish species inside marine protected areas (MPAs) and non‐MPA reference sites in the California Channel Islands, from 2003 (when MPAs were implemented) to 2017. We estimated harvest by fitting a size‐structured population model to survey data. Overall, harvest rates were effectively zero in MPAs but much higher in non‐MPA sites. This indicates successful adherence to MPA regulations, and possible displacement of fishing effort to reference sites. However, some poaching was detected in two MPA sites, highlighting the importance of assessing this quantity. This modeling approach could provide a tool to complement the long‐term management of MPA networks, particularly given the difficulty of acquiring harvest rate data at the spatial scale of individual MPAs.

of some types of fishing). Indeed, meta-analyses confirm that previously-fished populations do typically increase in abundance inside MPAs (e.g., Lester et al., 2009). However, a concern in any protected area is that continued, illegal harvest (i.e., poaching) will limit effectiveness (Ervin, 2003). Lack of enforcement of harvest bans can negate the positive effects of protected areas on target populations (Campbell et al., 2012;Sethi & Hilborn 2008;Watson, Dudley, Segan, & Hockings, 2014). For MPAs, Edgar et al. (2014) showed that successful enforcement of harvest regulations was one of the key features associated with successful protection.
Because poaching is-by definition-illegal, poachers have incentive to conceal their activities. Consequently, it is often difficult to quantify the intensity of poaching and its effect on population dynamics (Gavin, Solomon, & Blank, 2010). In marine fisheries, most data regarding poaching rates are qualitative rather than quantitative (e.g., anecdotal reporting rather than direct observation), hampering their utility (Bergseth, Russ, & Cinner, 2015). Further, the advent of global position system-based monitoring for fishing vessels has proven useful for spatial management of many commercial fisheries, but is difficult or impossible to implement in coastal, primarily recreational fisheries.
The potential for poaching inside MPAs is a problem for their successful adaptive management (note that in this paper we are exclusively referring to MPAs that prohibit harvest). Ideally, expectations for the magnitude and tempo of increase in fished populations are set based on life history factors and estimates of the pre-MPA level of fishing (Kaplan et al., 2019;White et al., 2011White et al., , 2013. Those expectations can then be compared to monitoring data to evaluate success (Grafton & Kompas, 2005;Nickols et al., 2019). If populations fall short of expectations, managers must consider the possibility of inappropriate MPA design (e.g., too small), mis-specification of the structure or parameters of the model used to set expectations, high variability in monitoring data, or noncompliance by fishers. Having an independent quantification of poaching improves that adaptive management decision process. Additionally, the possibility of poaching complicates other management efforts that rely on the assumption that MPAs are truly no-take (e.g., methods to estimate harvest in fished areas using MPAs as a reference point; Willis & Millar, 2004); directly testing that no-take assumption is thus critical.
Here, we describe a new ecological approach to quantifying poaching inside no-take MPAs. One effect of harvest on most marine populations is a truncation of the size structure: there are fewer larger (and older) individuals. Thus one can estimate harvest rates from the observed size structure of a population. White et al. (2016) developed a modeling tool to make such estimates from time series of size-structure observations. We applied that method to data from MPAs implemented in 2003 in southern California, USA (Figure 1), and found that overall, postimplementation harvest rates of four fish species were very low inside MPAs and quite high at non-MPA reference sites, indicative of high compliance with MPA regulations and potential redistribution of harvesting effort.

California MPAs and the Channel Islands marine reserve network
A network of 11 nearshore MPAs was established at the northern Santa Barbara Channel Islands in 2003 ( Figure 1; Botsford, White, Carr, & Caselle, 2014). The 11 MPAs comprise a mix of no-take reserves and conservation areas in which take of some species is permitted, though none permit any take of the species we considered in this study.

Visual diver survey data
Fish assemblages were surveyed annually from 1999 to 2017 at nearshore subtidal rocky reef sites throughout the northern Channel Islands by the Partnership for Interdisciplinary Study of Coastal Oceans (Caselle, Rassweiler, Hamilton, & Warner, 2015; Figure 1). Details of underwater visual survey protocols are in the Supplemental Information. The sites used in this study included locations inside MPAs (MPA sites) and locations in "reference" sites where harvest was allowed ( Figure 1; Table S1).

Study species
We focused on four fish species common on subtidal rocky reefs in the study area: blue rockfish (Sebastes mystinus, Scorpaenidae), kelp rockfish (Sebastes atrovirens, Scorpaenidae), California sheephead (Semicossyphus pulcher, Labridae), and kelp bass (Paralabrax clathratus, Serranidae). All four are targeted by recreational fishing and there are commercial fisheries for blue rockfish and California sheephead. These species have a relatively small home range size as adults, but a dispersive pelagic larval stage, so there is likely a high degree of demographic connectivity among MPA and non-MPA sites in the region (Mitarai, Siegel, Watson, Dong, & McWilliams, 2009

State-space integral projection model (IPM)
We represented fish population dynamics using a statespace IPM, a size-structured approach that allows for a direct comparison to size-abundance data. We followed the methods of White et al. (2016) to fit the model to the time series of size-abundance data for each species in each site, and estimate the harvest rate in the site based on the expected size distribution in an unfished population. Details of the model construction and analysis are provided in the Supplemental Information. The main parameter estimated by the model for each species-site combination was F, the harvest rate from 2003-2017 (after MPA implementation); this represents the harvest rate in non-MPA sites and the poaching rate in MPA sites. F is an instantaneous rate (units are y −1 ; the proportion of fish that survive harvest annually is e -F ); an easy shorthand is that F is approximately the proportion of the population harvested each year (as long as F is not much greater than ∼0.3).
We analyzed species-site combinations that had sufficient data for reliable estimation. Parameter estimation was done in a Bayesian manner, so we report posterior esti-mates of the harvest rates. Prior estimates for the harvest rates were based on published literature or stock assessments (Table S2). All model runs were inspected for convergence and checked for posterior predictive diagnostics (Table S3).

RESULTS
There were sufficient data for harvest rate estimation and models passed posterior diagnostics checks at two to eight individual study sites (both MPA and non-MPA) for each species (Table S1, S3; Figure 2). The two rockfish species are found in greater numbers in the cooler waters of the western islands; California sheephead and kelp bass are more abundant near the warmer-water eastern islands, hence the geographic pattern in data availability (Table S1; Figure 2). Posterior estimates of the post-2003 harvest rates (F, units y −1 ) were nearly zero for all four species in all of the MPA sites ( Figure 2). There were two exceptions: kelp rockfish and blue rockfish at Gull Isle and Painted Cave, respectively (both on Santa Cruz Island), had an estimated F that was greater than zero but less than 0.1 y −1 . By F I G U R E 2 Posterior distributions of harvest rate estimates (F) at each study site for (a) kelp rockfish, (b) blue rockfish, (c) California sheephead, and (d) kelp bass. Results for each site are displayed on a separate row, colored to denote protection status: red = MPA, blue = non-MPA. Sites with insufficient data for a given species are left blank. On each colored row, the posterior distribution of F is displayed as a gray histogram, with the median value indicated by a vertical black stripe. Note that for many of the near-zero estimates, the posterior distribution was very narrow and is obscured by the black median stripe. The final row in each panel displays the prior distribution on F (gray histogram). The inset in (c) shows the result for California sheephead at the Valley site, which required a different scale of the horizontal axis contrast, posterior estimates of the harvest rate at reference sites on Anacapa Island and the northern and southeastern shore of Santa Cruz Island were much greater than zero for multiple species (Figure 2). In most cases, the posterior estimates of the harvest rates at those sites were up to 10 times greater than the prior estimates of the harvest rates. By contrast, harvest rate estimates were near-zero at reference sites on the western shore of Santa Cruz island (Forney) and on the two western islands, Santa Rosa, and San Miguel.
Some of the estimated harvest rates were greater than one would typically expect to be possible, such as a posterior mean value of F = 2.7 y −1 for California sheephead at Valley (Santa Cruz Island). For reference, the prior estimate of the harvest rate for that species, based on a dataintegrated structured population model, was F = 0.219 y −1 for the Channel Islands region and F = 0.082 y −1 for the Santa Cruz Island population (Hamilton et al., 2011); a value of F greater than 1 implies that nearly all legal-sized fish at that site are harvested annually. Nonetheless, the F I G U R E 3 Examples of the model fit for kelp bass at two sites on Santa Cruz Island; (a) Cavern Point (inside an MPA) and (b) Pelican (not in an MPA). Panels show the model fit (black curve) to size abundance data (gray bars) for selected sample years near the beginning of data collection (2006)(2007) and at the end of the analyzed time series (2015)(2016)(2017). Upper and lower dashed lines indicate the middle 50% of the posterior predictive distribution of the size distribution, obtained by resampling the joint posterior parameter distributions. Population density is expressed as the abundance of fish per 30 × 2 × 2 m survey transect. Note that the vertical axis scale varies across panels to accommodate shifts in overall abundance and recruitment pulses. Note that the data (gray bars) are discrete counts in particular size bins, while the integral projection model output (black curves) represents a continuous population density distribution over size. Consequently, the bars and curve are sometimes different heights, but both model and data have similar integrals, that is, total abundance model using that posterior mean produced a good fit to the data, capturing interannual shifts in recruitment and adult size structure (Figure 3a). We inspected model outputs for all sites and species to ensure that they represented biologically plausible fits to the observed size structures (Figure 3b; Figure S1). Posterior predictive checks also confirmed that the models could recreate properties of the observed population dynamics (in this case, the average size of fish in the population each year) for most, but not all, species-site-year combinations ( Figure S2; note that sites with extreme deviations were excluded from this analysis and not shown in Figure S2). In some site-year combinations the posterior predictive model over-or under-predicted the average sizes.

DISCUSSION
Our analysis revealed that harvest rates are effectively zero inside most Channel Islands MPAs for four species frequently targeted by recreational and/or commercial fishers. By contrast, we estimated remarkably high harvest rates-up to 10 times greater than previous regional-scale estimates-at several sites where fishing remains legal. To our knowledge, this is the first study to directly assess whether MPAs have eliminated harvest in a fished population, rather than determining whether fish abundance or biomass increased (which could still be possible even if some poaching was occurring) or using surveys of fisher abundance in or near MPAs (e.g., Zellmer, Burdick, Medel, Pondella, & Ford, 2018a). This is an important advance in MPA assessment, because a major concern in MPA design is that fish will tend to move across MPA boundaries and be captured, even if the center of their home range is inside an MPA (Moffitt, Botsford, Kaplan, & O'Farrell, 2009;Zeller & Russ, 1998). That movement would reduce the effectiveness of the MPA, although in some cases such spillover is a desirable outcome to support fisheries. Our analysis suggests that such spillover is minimal for most of the Channel Island MPAs, because the fish populations within MPAs resemble unfished populations, not populations in which large individuals are removed if they stray into fished areas.
The spatial variation we detected in harvest rates at the reference sites aligns with our expectations about fisher behavior: fishing effort appears to be concentrated at Anacapa Island and the northeastern shore of Santa Cruz island, the sites closest to the ports where most fishing boats in the area are based (Figure 1), and sea conditions are typically more favorable for small boats in the eastern channel than further west (Zellmer, Claisse, Williams, & Pondella, 2018b). That we only found high harvest rates at the sites likely to be frequented by fishers lends additional support for the reliability of our estimates of zero fishing inside the MPAs on those islands.
The two exceptions to the overall pattern of nearzero harvest inside MPAs were low but non-zero harvest rates for kelp rockfish at Harris Pt. Reserve MPA on San Miguel Island and blue rockfish at Painted Cave MPA, on the northeastern shore of Santa Cruz Island. This result implies either that there is poaching in those MPAs, or that fish at those locations also stray outside of the MPA boundary. No poaching was detected for kelp bass or California sheephead at Painted Cave, but no other species was estimated at Harris Pt. Reserve. Based on that evidence we cautiously suggest the straying/spillover explanation is the most likely. Nonetheless, these results can guide managers in allocating enforcement resources toward particular locations. Indeed we recommend ongoing monitoring of activity at those sites.
One caveat in the interpretation of our results is that the shape of population size distributions depends on the interaction of growth and mortality rates, so size-based estimates of fishing mortality are sensitive to the values assumed for natural mortality and growth . We used empirically-based estimates of those parameters for all four species, but there is potential for spatial and temporal variation in both growth and mortality among our study sites, or between our study sites and populations and those used in the original literature estimates we referenced (Caselle et al., 2011). This could be particularly acute given the spatiotemporal variation in ocean temperature and other environmental variables among our study sites (Claisse et al., 2018;Hamil-ton, Caselle, Malone, & Carr, 2010). As such it would be worthwhile to re-evaluate our analyses if finer-scale estimates of fish growth and mortality rates became available. One challenging aspect of assessing the dynamics of species in MPAs is that reliable estimates of mortality and growth rates require killing large numbers of fish to extract otoliths for aging, which is an unappealing and selfdefeating proposition.
When interpreting the fits of model to data and the posterior predictive evaluations (as in Figure 3, S2), there were some instances in which the model predictions diverged from observed quantities in some years. This could lead to concern about model mis-specification. However given the nature of the data we believe a more likely explanation is extreme measurement variability in the underlying data, particularly because there was not a consistent trend of over-or under-prediction bias. Note that the model sometimes predicts the presence of small fish (young-ofthe-year) that were not observed in the data (these are often the most difficult size class to count), but appear in the data as the next-larger size class in the following year (e.g., for kelp bass at Pelican in 2016; Figure 3b). This phenomenon also occurred in a less-extreme manner for large size classes. This likely reflects the extreme variability inherent in annual underwater surveys: it is impossible for a single 1-2 day annual survey to produce a perfect sample of the full size distribution, particulary given variable physical conditions underwater. This reflects the need for a state-space approach applied across multiple years of data to account for missed observations (effectively, measurement error).
The estimated harvest rates outside of the MPAs in our study were considerably higher than prior estimates based on the most recent assessments of those stocks. This could be interpreted as evidence for fishing effort having been displaced from MPAs and concentrated in the remaining habitats, as is predicted to occur (Halpern, Gaines, & Warner, 2004) and has been observed in other MPA networks (e.g., Horta e Costa et al., 2013). We were not able to directly test that hypothesis because we had at most only four years of site-specific survey data prior to 2003 at these sites, and our simulation analyses have shown that amount of data is inadequate to obtain an unbiased estimate of harvest rates. However, the pattern of higher fishing at those sites is consistent with an observed increase (via aerial monitoring) in the number of recreational fishing vessels on the north shore of Santa Cruz Island and the south shores of Anacapa Island in the 2003-2008 period (Cabral, Gaines, Johnson, Bell, & White, 2017;Zellmer et al., 2018a). The higher estimated fishing rates for kelp bass relative to the other three species are also consistent with kelp bass being the top species caught by recreational fishers in Southern California since the 1930s. What remains to be investigated is the degree to which fishery yields might increase as a result of larval spillover from MPAs, if effort in fished areas has increased accordingly (Halpern et al., 2004;Kerwath, Winker, Götz, & Attwood, 2013).
The "success" of any given MPA depends on the stated goals of that MPA. In the case of the Channel Islands MPAs, the primary goal was to conserve ecosystem biodiversity (Botsford et al., 2014). Indeed, recent studies have found an overall increase in biomass of fishery-targeted species in Channel Islands MPAs (Caselle et al., 2015). However, there are limitations to using organism abundance data as a sole indicator of MPA management success. First, it is possible to observe increases in abundance in MPAs even if poaching is reduced but not eliminated. Second, interannual variability in larval recruitment can delay increases or make abundance trends highly variable and difficult to detect. Nickols et al. (2019) recently showed how both factors can mask eventual positive trends in fished populations in MPAs in central California. Additionally, increases in fish biomass in MPAs are often assessed relative to biomass in fished areas, but it is possible for that inside-outside ratio to be large regardless of whether the population in the MPA is increasing or decreasing (Moffitt, White, & Botsford, 2013), making it difficult to assess MPA performance using "control" sites. Our approach of using size distributions to estimate harvest rates directly eliminates the confounding effect of recruitment variability, and allows us to confirm that MPA regulatory enforcement is actually successful.
Our results will inform the pending assessment of California's MPA network in 2022 (California Department of Fish and Wildlife & California Ocean Protection Council, 2018) and offer a potential solution for estimating poaching rates inside MPAs and harvest rates outside of MPAs. This type of information is lacking for most MPA locations around the world, but is critical for assessing effectiveness and for adaptive management (Bergseth et al., 2015;Brown et al., 2018). Previous works have clearly shown a relationship between compliance and MPA success (Aburto-Oropeza et al., 2011;Giakoumi et al., 2017). Yet many MPA managers are resource-limited and unable to monitor human uses in and around MPAs, even in cases where biological monitoring exists. Importantly, our model does not require extensive, professionally-collected data on fish or fishers, as the inputs are simply a time series of fish numbers and sizes. With the increasing availability of inexpensive, video cameras configured in stereo, fish sizes can be readily and accurately surveyed by stakeholders or citizen scientists (Lowry, Folpp, Gregson, & Suthers, 2012), or potentially using less labor-intensive remote methods, such as remote operated vehicles (Perkins et al., 2020) or baited remote underwater video (Langlois et al., 2020). Elsewhere we have investigated the data requirements for our method, and discovered that unbiased estimates of harvest rates require observations of approximately 100 individuals per year for seven years (with fewer years possible if the size distribution is better characterized; Yamane, L., et al., unpublished manuscript). This level of sampling could be daunting for some monitoring programs (and it would be preferable to have information before seven years had elapsed). Fortunately, that analysis also indicates that even with much poorer datasets (∼50 individuals per year for three to four years), one could detect a poaching rate of at least F = 0.05 y −1 , so inferences would still be possible on a coarser scale.
This approach provides a novel alternative to expensive human-usage surveys such as drones or aerial surveys, and can be implemented at scales appropriate to small MPAs, unlike the majority of fisheries landings data. Poaching is a global threat to the effectiveness of MPAs (Giakoumi et al., 2018), so length-based methods leveraging ecological monitoring data such as this one offer a way forward. The method we have developed could also be applied to other harvested species with indeterminate growth in which size distributions are an effective indicator of harvest, such as freshwater or marine turtles (e.g., Norris, Peres, Michalski, & Gibbs, 2019;Piacenza, Balazs, Hargrove, Richards, & Heppell, 2016).
Our assessment of poaching rates in MPAs was only possible because of the availability of long-term surveys of fish populations at these islands. A key component to adaptive management is "learning", which is difficult (perhaps impossible) without long-term monitoring of the populations that respond to that management. There are some creative options based on fishery-oriented approaches (e.g, Harford et al., 2016), but in general it is essential to continuously monitor protected areas-both marine and terrestrial/freshwater-if one desires to ensure they are achieving their stated goals (Danielsen et al., 2005;Grorud-Colvert et al., 2014;White et al., 2011). Such long-term data can not only contribute to an understanding of the ecosystems protected by MPAs, but as demonstrated here, can also provide much needed information on harvest rates and compliance.

A C K N O W L E D G M E N T S
This research was supported by the National Science Foundation (NSF: OCE-1909303, OCE-1820540, OCE-1758000). We thank the legion of divers who collected the data used in this analysis. This is publication 513 of the Partnership of Interdisciplinary Study of Coastal Oceans, funded primarily by the David and Lucile Packard Foundation (2019-68104).