Collaborative fisheries research reveals reserve size and age determine efficacy across a network of marine protected areas

A variety of criteria may influence the efficacy of networks of marine protected areas (MPA) designed to enhance biodiversity conservation and provide fisheries benefits. Meta‐analyses have evaluated the influence of MPA attributes on abundance, biomass, and size structure of harvested species, reporting that MPA size, age, depth, and connectivity influence the strength of MPA responses. However, few empirical MPA evaluation studies have used consistent sampling methodology across multiple MPAs and years. Our collaborative fisheries research program systematically sampled 12 no‐take or highly protective limited‐take MPAs and paired fished reference areas across a network spanning 1100 km of coastline to evaluate the factors driving MPA efficacy across a large geographic region. We found that increased size and age consistently contributed to increased fish catch, biomass, and positive species responses inside MPAs, while accounting for factors such as latitude, primary productivity, and distance to the nearest MPA. Our study provides a model framework to collaboratively engage diverse stakeholders in fisheries research and provide high‐quality data to assess the success of conservation strategies.


INTRODUCTION
Conservation strategies such as marine protected areas (MPAs), including no-take marine reserves, are popular spatial management tools used to protect marine biodiversity, promote resilient ecosystems, and deliver benefits to society (Mouillot et al., 2020).Due to the success of many MPAs that limit fishing activities in supporting species abundances (Cinner et al., 2018) and human well-being (Ban et al., 2019), policymakers have also leaned on these strategies to help mitigate the effects of global climate change (Naeem et al., 2022).However, the measured outcomes of MPAs can vary greatly across space and time based on protection level, enforcement, MPA age, and various ecological conditions (Gill et al., 2017;Grorud-Colvert et al., 2021;Meehan et al., 2020).Global meta-analyses suggest that the cessation of fishing in MPAs generally increases the abundance of targeted fish species (Edgar et al., 2014;Giakoumi et al., 2017), with more delayed indirect effects on nontargeted species (Babcock et al., 2010).However, for any single MPA, the abundance of fished species may remain constant or even decrease through time, possibly in response to external factors, such as oceanographic conditions, recruitment variability, changes in food webs or habitats (Edgar & Barrett, 2012) and gravity of human impacts (Cinner et al., 2018).Such variable results have led researchers to investigate the specific factors that drive organismal responses in MPAs (Nickols et al., 2019) along with increased development of MPA networks, which are designed to buffer populations from those external drivers, rather than stand-alone MPAs (Grorud-Colvert et al., 2014).
Theoretical benefits of MPA networks include protecting diversity of habitats and species, reducing localized impacts of climate change, maintaining genetic diversity and population persistence (White et al., 2011), and creating a more equitable distribution of fisheries costs and benefits (Leenhardt et al., 2015).Understanding how species respond to MPA protection is essential to set appropriate expectations for stakeholders, as well as to adaptively manage the network into the future.Theoretical studies and meta-analyses suggest that reserve characteristics such as size, age, connectivity, location, and adjacent fishing effort may influence the magnitude of positive effects from MPA networks (Claudet et al., 2008;Edgar et al., 2014;Lester & Halpern, 2008).While meta-analyses can provide a wealth of information on diverse ecosystems and species across space and time, the data quality can vary greatly (Halpern & Warner, 2003) and empirical data are often collected in using different methods (e.g., visual census, fishing, underwater video; Claudet et al., 2010), influencing the conclusions researchers draw.
From 2003 to 2012, the government of California, a state located on the Pacific coast of the USA, implemented a statewide network of MPAs to restrict harvest of marine species and to provide conservation and fisheries benefits using the best available science on habitat, species diversity, dispersal capability, and distribution of human uses.This process resulted in a scientifically designed MPA network that includes 124 MPAs, spanning 1100 km of coastline and 2200 km 2 (16%) of state waters (Botsford et al., 2014).Baseline monitoring of the newly established network was conducted across diverse habitats, including kelp forests and rocky reefs.Data from these monitoring programs indicated that both MPAs and areas open to fishing within the same geographic area had similar abundance, biomass, and species composition when MPAs were first enacted (e.g., North coast: Mulligan et al., 2017;Central coast: Starr et al., 2015, South coast: Caselle et al., 2015).
We evaluated the effects of no-take or limited take reserves at 12 MPAs within the California MPA network using 6 years (2017-2022) of standardized hook and line surveys conducted by the California Collaborative Fisheries Research Program (CCFRP).CCFRP is a unique partnership among scientists, managers, and anglers that combines the expertise and ideas of diverse stakeholders to gather information on fish stocks for management and conservation (Wendt & Starr, 2009) while improving buy-in for fisheries management (Mason et al., 2020).By using catch data collected by CCFRP inside and outside MPAs, we are able to assess how a variety of metrics including MPA age, size, adjacent fishing pressure, geographic location, network connectivity, and environmental factors affect total catch rates, catch biomass, and the number of species exhibiting positive responses to MPA protection (i.e., higher biomass inside MPAs compared with reference sites).We expected catch rates and catch biomass to be higher in larger and older MPAs (McClanahan & Graham, 2015).Further, we expected more rapid population responses and stronger MPA effects at lower latitudes in southern California where human population density (Free et al., 2023) and fishing pressure are greater and warmer water temperatures are associated with faster growth rates and more consistent recruitment of the dominant taxa (Caselle et al., 2023).

CCFRP surveys
CCFRP collaborates with volunteer anglers and the fishing industry to conduct surveys of rocky reef-associated fish communities inside and outside of MPAs using standardized, repeatable hook-and-line fishing protocols (Wendt & Starr, 2009) S3).Each of the 24 total sites was sampled with standardized hook-and-line protocols for 3 days (i.e., visits) per year, with four randomly selected grid cells (each 500 m × 500 m in area) sampled per day (Starr et al., 2015).Within each grid cell, we completed three 15-min fishing periods, where volunteer anglers actively fished over rocky reef habitat.All fish captured were identified to the lowest taxonomic level, measured to total length (cm), and released.For additional details on CCFRP, reference site selection, and fishing protocols please refer to the Supporting Information Text 1.

MPA characteristics
We calculated the age of each MPA by subtracting the year of implementation from the sampling year.MPAs were implemented in various years: Channel Islands in 2003, central coast in 2007 and northern and southern (non-island) coasts from 2010 to 2012.Total protected area (km 2 ) was defined as the MPA area for solitary MPAs and the sum of protected areas for paired MPAs (locations where no-take SMRs and highly protective limited-take SMCAs were adjacent to each other).We extracted the latitude for the centroid of each MPA for spatial analyses.
We calculated the distance from the nearest MPA (excluding SMR/SCMA pairs) within the network to assess the influence of connectivity and nearest port as a proxy for fishing pressure (Harborne et al., 2018).These distances were measured as the shortest linear distance over water from the boundary of the MPA to the nearest MPA or port using ArcGIS PRO (version 3.0.3).Across our study system, latitude and human population density are highly correlated and previous work along this coast suggests that human population density may not be the best indicator of human engagement with MPAs (Free et al., 2023), so we did not consider that metric in our MPA characteristics.Fishing pressure directly adjacent to an MPA has been shown to modulate the magnitude of positive effects on fish populations (biomass and catch) across a portion of this network (Ziegler et al., 2022), but recreational fishing data at the statewide-scale are limited and less reliable in some regions.Therefore, we calculated an index of recreational fishing pressure outside each MPA as the total number of groundfish landed between 2017 and 2019 for the nearest port complex (from the CalFISH database; Free et al., 2022) divided by the distance from the nearest port (under the assumption that realized fishing pressure declines with distance from port).Due to lack of highquality habitat mapping across all of our sampling sites, we used the habitat rugosity metric from each CCFRP fishing drift (Supporting Information

Analysis
We calculated both catch per unit effort (CPUE) and biomass per unit effort (BPUE) as relative metrics of abundance and biomass in each MPA and REF site.CPUE was calculated by dividing the total number of fishes caught by total angler hours of fishing in all grid cells in a day (catch angler h −1 ).For each individual fish, we converted length to biomass using published length-weight relationships (Love et al., 1990;Froese & Pauly, 2023).We then calculated BPUE as the total weight of fish caught in kilograms divided by the number of angler hours fished (kg angler h −1 ).CPUE and BPUE for the grid cells sampled on a given day was then averaged to estimate CPUE or BPUE for the total fish community inside and outside each MPA for a given sampling year.Data were tested for normality with a Shapiro-Wilk test (W = 0.94, p < 0.001).
To assess the drivers of spatial variation in the strength of MPA responses across the network, we calculated a yearly response ratio for each MPA using both CPUE and BPUE values (Stewart, 2010).Response ratios were calculated as the log of the value (X) inside the MPA divided by the value (X) at the REF for year i: ) A ratio above 0 indicates higher CPUE or BPUE inside the MPA relative to REF sites, while the magnitude of the ratio reflects the relative difference in CPUE or BPUE between each MPA and its associated REF (i.e., the strength of MPA response) on the log 10 scale.Using a BPUE response ratio calculated independently for each species, we determined the number of positive species responses for each MPA and year.
To determine the relative influence of MPA characteristics on the efficacy of MPAs to enhance fish catches (i.e., response ratios), we fit generalized linear mixed effects models with normally distributed errors using the lme4 package in R (Bates et al., 2015; Table 1).Prior to anal-yses, we tested for collinearity among fixed effects using a Pearson correlation matrix (Figure S4).Both SST and the difference in relief between MPAs and REFs were highly correlated with latitude (r > 0.7) and therefore were removed from subsequent models.Our final models included MPA age, latitude, total protected area (km 2 ), relative fishing pressure (landings [n] per distance from nearest port [km]), mean NPP (mg C m −2 day −1 ), and distance to the nearest MPA (km) within the MPA network as fixed effects with a random effect of sampling site.We also considered interactions that were deemed the most ecologically relevant or plausible: MPA age and latitude, fishing pressure and latitude, and total protected area and fishing pressure.For the metric of the number of positive species responses, including the random effect overfit the models and therefore we used negative binomial models with no random effects.Akaike information criterion corrected for small sample sizes (AICc) was used to compare models with all possible variable combinations.Models were ranked using their AICc difference from the model with the minimum AICc (ΔAICc), and we considered all models with ΔAICc < 2 for interpretation.Beta coefficients (β) were standardized by standard deviation for direct comparisons across models.We then conducted model averaging in the MuMIn package (Bartoń, 2023) to assess which variables best predicted the response variables across the top models and assessed the partial effects of each factor on MPA responses (Burnham & Anderson, 2002).All analyses were conducted in R statistical software version 4.3.0(R Development Core Team 2023).

RESULTS
We observed greater total catch rates (CPUE) and catch biomass (BPUE) inside MPAs compared with REFs at all locations surveyed across the network of MPAs from Northern to Southern California (Figures 1 and S5, S6).Both CPUE and BPUE response ratios varied across the MPA network through time (Figure 2a).Of the 12 MPAs sampled, 7 had positive CPUE response ratios and 10 had positive BPUE response ratios in all sampling years (Figure 2b), indicating consistently higher total catch rates and greater biomass inside MPAs.The largest effect sizes for both CPUE and BPUE (100-fold higher BPUE in MPA relative to the REF) occurred at Swami's SMCA (where absolute values were low in both MPA and REF) and at Anacapa Island SMR/SMCA (where values for both CPUE and BPUE in REF were extremely low) (Figure 2).Overall, 19 of 72 possible year-site combinations for CPUE and 29 of the combinations for BPUE showed response ratios of >=1, respectively (10-fold higher values of catch and biomass inside MPAs relative to REFs).
TA B L E 1 MPA characteristics for each of the 12 MPAs sampled by CCFRP.MPA characteristics include year of MPA establishment, MPA name, whether the MPA is solitary (no) or adjacent to another marine managed area with less restrictive no-take regulations (yes), the amount of total protected area (km 2 ), the latitude of the centroid for the MPA, the nearest port name, distance to the nearest port in kilometers, groundfish landings (no.groundfish y −1 ) extracted from the CalFish database (Free et al., 2022), fishing pressure (landings/distance to nearest port), distance to the nearest MPA in kilometers and the difference in relief calculated by dividing the mean rugosity score (1-3) from each fishing drift in a MPA by the rugosity score (1-3) from each drift in the REF.Mean net primary production (mg C m −2 day −1 ) and sea surface temperature ( • C) are not shown as they were examined on an annual basis.Four competing models best fit the CPUE data, including variables of MPA age, area, latitude, mean NPP, and distance to the nearest MPA.After model averaging the predictor variables that best explained CPUE response ratios were MPA age (with coefficient, β = 0.29, p = 0.006) and total protected area (β = 0.50, p = 0.01; Figure 3).For BPUE response ratios, six competing models best fit the data, including variables of MPA age, area, latitude, mean NPP, distance to the nearest MPA, and fishing pressure.However, similar to CPUE, MPA age (β = 0.29, p = 0.002) and total protected area (β = 0.29, p = 0.002) were the top predictors of BPUE response ratios across all MPAs.MPA age and total protected area were positively associated with CPUE and BPUE responses ratios observed (Table 2), such that MPAs that were both larger and older tended to have higher total catch rates and fish biomass inside MPAs relative to the reference sites (Figure 4).

Region
The number of positive species responses varied across MPAs and as a function of sampling year, with the percentage of positive species responses ranging from 40 to 100%.
Of the 72 possible year-site combinations, 65 had more positive species responses than negative species responses (>50%; Figure 5).Three models best explained, the total number of positive species responses within MPAs over time.After model averaging, the top predictor variables were total protected area (β = 0.12, p < 0.001), fishing pressure (β = 0.45, p = 0.03), as well as interactive effects of area and fishing pressure (β = −0.50,p < 0.001; Table 2).The number of species exhibiting positive responses to MPAs was higher in larger MPAs and in locations where fishing pressure was higher outside the MPA; however, the fishing pressure effect was dampened in locations with large MPAs.

DISCUSSION
Size and age of no-take or highly protective MPAs positively affected fish catch rates and catch biomass.For many of the MPAs sampled, response ratios were greater in 2022  compared with 2017 when CCFRP expanded sampling across the state (the program started in central California in 2007; Starr et al., 2015), even though the MPAs were implemented at different times in different regions between 2003 and 2012.Extensive efforts (such as consulting multibeam habitat maps and members of the charter fishing community) were made at the beginning of the CCFRP program to identify references areas that were equivalent in habitat to the MPAs sampled (Yochum et al., 2011).Similarly, baseline sampling and times series analysis is select regions indicated that fish density and biomass values were similar among MPAs and REFs at the time of MPA implementation across the California coast (Caselle et al., 2015;Hamilton et al., 2010;Starr et al., 2015;Mulligan et al., 2017;Ziegler et al., 2022).As a result, we interpret that the differences we observed in catch rates and biomass are due to the removal of fishing pressure inside the MPAs and not pre-existing differences.
In general, we observed stronger responses to MPA protection in older MPAs (Anacapa, Carrington Pt., Pt.Lobos) compared with younger MPAs implemented in 2010-2012, closer to the expansion of CCFRP to a statewide MPA monitoring program.Synthetic analyses of European MPAs similarly show that fish density increases with the age of the MPA (Claudet et al., 2008) while global analyses report that older MPAs (>10 years since implementation) host increased fish biomass (Edgar et al., 2014).However, we did observe comparably high responses in some newer MPAs in each region (e.g., Swami's, Stewarts Pt.), although these MPAs tended to have higher total amounts of protected area and Swami's had much lower BPUE overall.Many of the response ratios are relatively flat over the time series, possibly indicating that the effects of the MPAs have begun to saturate (McClanahan et al., 2007;Russ & Alcala, 2010), that extensive spillover may be occurring (Russ & Alcala, 2011) or that recovery is occurring very F I G U R E 5 Tornado plots of species-specific average BPUE response ratios for all MPAs surveyed between 2017 and 2022.Values > 0 indicate a positive species response (i.e., higher biomass inside MPA relative to REF), while values < 0 indicate a negative species response.Paired MPA/REF sites are ordered and colored by latitude (and relative differences in sea surface temperature) (blue = northern and cool, red = southern and warm).Bars are the mean response ratio per species and errors are ±95% CI.
slowly.McClanahan et al. (2007) reported that densities of coral reef fish in Kenya may reach an asymptote anywhere from 5 to 10 years or to up to 37 years after MPA protection, while Russ and Alcala (2010) demonstrated that predatory fish biomass saturated after 20-40 years of protection in the Philippines.While it is possible that total fish abundance and biomass in MPAs across this network may saturate relatively quickly, the studies that found rapid recovery examined tropical species that are likely to have much higher growth rates than those in temperate waters in California.
We also found that the differences in fish CPUE and BPUE between MPAs and REFs were greater in locations where the MPAs were larger.Theoretical studies have hypothesized that larger MPAs would be more effective at increasing biodiversity (Krueck et al., 2018) and density of economically important species (Krueck et al., 2017).Past empirical studies have provided evidence that larger MPAs are characterized by greater increases in fish density and biomass (Claudet et al., 2008;Edgar et al., 2014).Larger MPAs may allow a higher proportion of fishes with larger home ranges to remain protected within the MPA, compared with smaller MPAs, where more mobile species may be likely to spillover into unprotected areas and be caught by anglers (Di Lorenzo et al., 2020).In contrast, other largescale global syntheses of MPAs are more equivocal on the importance of MPA size (Guidetti & Sala, 2007) or report greater effects in smaller MPAs (Giakoumi et al., 2017).
A recent meta-analysis reported that edge effects in MPAs can effectively reduce the MPA area (Ohayon et al., 2021).The authors concluded that fishing-the-line behavior-where fishers disproportionately fish the edge of a protected area-may be largely responsible.These MPA edge effects may be particularly pernicious for small MPAs (<10 km 2 ), since the total area to perimeter ratio is small, providing more opportunities for fishing-the-line relative to the size of the MPA.Larger MPAs may also increase larval retention, promoting self-recruitment and the build-up of biomass inside the MPA boundary, increasing overall MPA efficacy (Krueck et al., 2017).Our results indicated that catch biomass was higher inside every MPA compared with its associated REF, with MPA effectiveness increasing further with increased with MPA size.Thus, while small MPAs are still effective at protecting fish populations, our statewide data revealed that larger ones within the network were even more effective.
Overall, we expected more rapid population responses and stronger MPA effects with decreasing latitude in the warmer waters of Southern California (Caselle et al., 2023) due to faster growth rates and higher recruitment of the dominant southern taxa (kelp bass, California sheephead, ocean whitefish) compared with northern taxa that generally have slower growth and more sporadic recruitment (rockfishes and lingcod) (Caselle et al., 2010).For example, Hamilton et al. (2010) observed rapid increases in biomass of southern taxa following the establishment of the Channel Islands MPAs, while Starr et al. (2015) and Ziegler et al. (2022) reported much slower responses of rockfishes in the initial years following MPA establishment in Central California.Consistent with our hypotheses, MPA responses tended to decline with increasing latitude (ΔAICc < 2).Latitude is also negatively correlated with SST and our metric of fishing pressure, but positively correlated with NPP (Figure S4), making it challenging to disentangle which factor(s) are ultimately driving the latitudinal response to MPAs.
We also observed a weak positive effect of adjacent fishing pressure on the strength of MPA responses, especially the number of positive species responses; shown to be important in other California studies (Jaco & Steele, 2020;Ziegler et al., 2022) and globally (Griffiths et al., 2022).Fishing pressure is much more intense in southern California (∼1 million recreational anglers), and much of the coast is easily accessible.In contrast, human population density and recreational fishing pressure are lower in central and northern California and there are smaller, more disparate ports in these regions, creating more remote sections of the coastline that are much less accessible to boat-based anglers (Free et al., 2023).The limited strength of the relationship between fishing pressure and MPA response ratios in this study may be due to a variety of factors.Most importantly, the fishing pressure metric (e.g., landings / distance from port) was likely too coarse to provide strong predictive power at the individual MPA scale.Still, we observed a positive association between the strength of the MPA response and fishing effort, similar to Ziegler et al. (2022) on the central California coast, indicating that fishing pressure is positively correlated with MPA responses across the larger scale of the statewide network.Higher fishing pressure in references sites equates with bigger differences in catches and biomass inside MPAs compared with those fished areas.These results emphasize the need for high quality, spatially explicit fishing pressure data to better assess the effectiveness of MPAs across large spatial scales.
Overall, our results-using standardized data collected through a collaborative fisheries research program, CCFRP-support and expand much of the previous literature on the attributes that result in positive effects of MPAs to enhance abundance and biomass of targeted species.Our results show that targeted fish abundance increased in MPAs throughout the network, and that MPA size, age, and likely some combination of latitude, fishing pressure, and temperature had significant effects on the strength of the MPA response.Critically, this partnership of academic and agency scientists with the sportfishing industry and local anglers is able to generate high quality, reliable, and transparent information on the status of species, communities, and ecosystems protected by MPAs essential for adaptive management.While metaanalytical frameworks and studies can provide insight into global patterns of MPA efficacy, standardized sampling programs that engage stakeholders and span MPA networks can provide essential data for local and regional marine resource management.Our approach combining community engagement and a standardized and statistically rigorous sampling design not only provides valuable insight for the design and adaptive management of MPA networks spanning large geographic scales, it provides a model framework for the engagement of local stakeholders, often wary of conservation strategies, into MPA monitoring and fisheries management.

A U T H O R C O N T R I B U T I O N
a) The average catch per unit effort (CPUE; no angler h −1 ) and biomass per unit effort (BPUE; kg angler h −1 ) of fishes at 12 marine protected areas (MPAs, red bars) and reference sites (REFs, blue bars) along the coast of California across all 6 years sampled statewide (2017-2022).Values are mean ± 95% CI.(b) Map of 12 MPAs (red polygons) sampled along the California coast.
Catch per unit effort (CPUE; no angler h −1 ) and Biomass per unit effort (BPUE; kg angler h −1 ) response ratios for the 12 marine protected areas (MPAs) spanning 10 degrees of latitude along the coast California from 2017 to 2022.Paired MPA/REF sites in the legend are ordered and colored by latitude and relative differences in sea surface temperature (blue = northern and cool, red = southern and warm).Response ratio values (y-axis) greater than zero indicate higher total fish biomass inside the MPAs relative to the reference sites.Dashed line indicates no difference in CPUE or BPUE between MPA and REF sites.
Partial effects plot of the full model for CPUE response ratios with all variables included in model averaging.Partial residuals of response ratios as a function of (a) total amount of protected area (km 2 ), (b) MPA age, (c) latitude, (d) mean NPP, and (e) distance to nearest MPA while other fixed effects were held constant at median values.Trend lines indicate conditional fit of the models without the random effect of site.
Partial effects plot of the full model for BPUE response ratios with all variables included in model averaging.Partial residuals of response ratios as a function of (a) total amount of protected area (km 2 ), (b) MPA age, (c) latitude, (d) mean NPP, (e) distance to nearest MPA, and (f) fishing pressure, while other fixed effects were held constant at median values.Trend lines indicate conditional fit of the model without the random effect of site.
−2 day −1 ; NPP) from the Central and Northern California Ocean Observing System Repository.SST was collected by Advanced Very High-Resolution Radiometer instrument aboard NOAA's Polar Operational Environmental Satellites from 2004 to 2022 at a 1.47 km resolution and NPP was collected by the California Current Merged Satellite daily from 1996 to 2022 at a 4-km spatial resolution.The spatial average of both SST and NPP within the boundaries of each MPA were extracted annually for each year of CCFRP sampling.See Table1for a summary of MPA characteristics for all sampling sites.
by the average rugosity scores for drifts in the associated REF.Scores of 1 indicate the relief in the sampled areas of the MPAs and REFs were equal; while scores greater than 1 indicate higher relief in the sampled areas of the MPA compared with the REF.We also extracted environmental variables of sea surface temperature ( • C; SST) and net pri-mary production (mg C m Comparison of top eight models predicting the influence of MPA characteristics (MPA age, total protected area, relative fishing pressure, latitude, mean NPP and distance to the nearest MPA) on CPUE and BPUE response ratios, and number of positive species responses across 12 MPAs in a temperate marine MPA network.For each model, marginal R 2 TA B L E 2 c), degrees of freedom (df), log likelihood (LL), AICc, and ΔAICc, are shown.Shaded models had a ΔAICc < 2.