Diving deep into the network: Quantifying protection effects across California's marine protected area network using a remotely operated vehicle

Abstract

Networks of marine protected areas (MPAs) are being increasingly implemented worldwide as conservation management tools. We report here on MPA effectiveness using data from a 17‐year long remotely operated vehicle (ROV) monitoring program spanning mesophotic (~ 20–130 m) depths in 26 MPAs across California's MPA network. We utilize a spatial modeling approach that includes important environmental covariates as well as spatial dependence in the data, and allows the separation of statewide and regional trends in the abundance of focal species from additional trends specific to MPAs. We demonstrate that there have been statewide and/or regional recoveries in abundance for the majority of our 10 focal demersal fish species, with all four statewide species assessments and 18 out of 22 species‐region combinations assessed displaying positive trends. We also demonstrate that MPA protection has had an additional positive effect on the abundance of the majority of these focal species inside MPAs compared with reference areas, with positive effects for all four statewide species, and 18 out of 22 statewide/regions assessed showing positive effects, four showing no statistically detectable differences, and no negative MPA effects found. Comparisons with theoretical expectations of MPA recovery for our focal species showed that 2 out of 4 statewide, and 11 out of 17 species‐region combinations assessed displayed higher mean MPA effects than expected. Our results highlight strong trajectories of increasing abundance and additional MPA effects for many of our focal species, demonstrating that MPAs are having positive effects in mesophotic depths across the network as well as at previously reported shallower depths, and that image‐based platforms such as ROVs provide an important tool to support timely reporting on the effectiveness of MPA networks.

INTRODUCTION

Marine protected areas (MPAs) are an important conservation tool to mitigate human impacts on marine ecosystems (Edgar et al., 2014; Lester et al., 2009; Sala et al., 2021). Management goals often include the preservation of marine resources, the rebuilding of fisheries stocks, and the promotion of ecosystem integrity (Agardy, 2000; Grorud-Colvert et al., 2021; Klein et al., 2008). In many cases, these goals are mandated by legislation, which contains explicit conservation targets. Considerable importance is therefore placed on the timely reporting of MPA effectiveness in meeting stated goals. Metrics for quantifying MPA effectiveness often involve tracking trajectories of abundance and/or biomass of key-targeted species both inside and outside of MPAs. Reporting on MPA effectiveness, however, particularly across a large network of MPAs, requires a concerted effort to attain high-quality monitoring data.

Obtaining quality monitoring data in subtidal marine environments is typically expensive and technically challenging. Monitoring data sets are consequently sparse in space and time, making inference on temporal changes difficult due to the inherent uncertainty of limited data. Part of the technical challenge is that many important habitats for key species can span depths that lie beyond the limits of more traditional methods such as SCUBA surveys. Ecosystems in depths ~20–150 m, often referred to as mesophotic ecosystems (e.g., Cerrano et al., 2019), are being increasingly acknowledged as critical to our understanding of the dynamics of nearshore marine ecosystems as they often contain large areas of habitat for species of management interest (e.g, Castellan et al., 2022), such as previously fished species across an MPA network (Starr et al., 2022). Therefore, technical solutions are required that capture the dynamics across the depth ranges occupied by target organisms, and with sufficient spatial and temporal resolution to detect changes attributable to MPA establishment.

Image-based sampling platforms such as remotely operated vehicles (ROVs) and baited remote underwater video (BRUVs) provide a means of collecting nonextractive video data on the abundance and sizes of target organisms and the habitats in which they reside across large spatial scales and depth ranges (Knott et al., 2021; Sward et al., 2019). Evidence is beginning to mount that these platforms can be effective tools for detecting the effects of MPAs in mesophotic depths (e.g., Knott et al., 2021; Vigo et al., 2023). While ROVs have been used for some time in exploratory surveys of the sea floor, they are also being increasingly used in monitoring programs in deeper continental shelf waters. They are particularly suited to MPA monitoring as they provide a means of nonextractive sampling inside protected areas. To date, however, most studies using ROVs for MPA monitoring are limited in the temporal scope of reporting, with protection effects often inconclusive (e.g., Haggarty et al., 2016; Karpov et al., 2012). In this study, we utilize a 17-year long data set of ROV surveys conducted inside MPAs and nearby reference areas across mesophotic (~20–150 m) depths in the California MPA network with the explicit aim of quantifying the effect of MPA establishment. This data set provides a particularly timely opportunity to explore MPA effects in mesophotic depths across the California MPA network as (i) it is spatially extensive, covering a large number of MPAs, and (ii) these MPAs have now been in place between 11 and 17 years, which may be sufficient time to detect MPA effects for some key previously fished species in this study system (Kaplan et al., 2019; Perkins et al., 2020; Starr et al., 2015).

One of the expected effects of MPAs is an increased density of previously fished species as the length of protection increases (Claudet et al., 2008; Halpern & Warner, 2002; Hopf et al., 2016). The California MPA network provides an ideal study system for testing the effect of MPAs on the density of targeted fish species across large scales and depths due to the large number of MPAs that include mesophotic depths, and sufficient lengths of protection to detect these effects for previously fished species. Previous work across MPAs in California's Channel Islands has shown that MPAs are having positive effects for previously targeted fish species across shallower depths (< 20 m) in these longest established MPAs in the network (Caselle et al., 2015; Hamilton et al., 2010). However, a number of challenges exist, especially when assessing MPA effectiveness over large scales of management interest, such as California's MPA network. These include other historical and ongoing fisheries management measures, in addition to MPAs, that also impact on temporal trends. For example, the United States Pacific Fishery Management Council implemented a number of fisheries depth closures along the west coast in 2002 to help rebuild stocks of overfished species (Federal Register, 2003). Stock assessments on the west coast of the US have indicated trends for some species may be increasing in recent years (e.g., Dick et al., 2021; Fisheries, 2024; Johnson et al., 2021; Monk et al., 2021; Taylor et al., 2021; Wetzel et al., 2021a; Wetzel et al., 2021b), perhaps in part due to these spatial closures (e.g., Keller et al., 2019) or other fisheries management measures such as commercial and recreational catch limits. For fisheries and MPA management purposes, there is therefore the need to be able to separate regional and statewide trends in abundance that may be due to other fisheries management measures from what may be occurring inside MPA boundaries.

The focus of stock assessments and MPA monitoring is often individual species as they differ in many ways, including in their life history characteristics, with consequent levels of susceptibility to impacts resulting in different trajectories of recovery or decline. In the case of California's MPA network, previous work has highlighted that estimates of fishing mortality prior to MPA establishment combined with a knowledge of the life history traits of individual species can allow projections of expected magnitudes and timelines of recovery (Kaplan et al., 2019). These expected trajectories of recovery for individual species can be compared with monitoring data to ascertain whether MPAs are meeting expectations (Nickols et al., 2019).

Here, we use a spatially and temporally extensive data set of ROV surveys to test the effect of the length of MPA protection on the density of 10 focal species that are key-targeted species by recreational and commercial fisheries: brown rockfish (Sebastes auriculatus), California sheephead (Semicossyphus pulcher), copper rockfish (Sebastes carinus), gopher rockfish (Sebastes carnatus), kelp greenling (Hexagrammos decagrammus), lingcod (Ophiodon elongatus), quillback rockfish (Sebastes maliger), vermilion rockfish (Sebastes miniatus), and yelloweye rockfish (Sebastes rubberrimus). These species were chosen as they are all demersal groundfish species that are captured well by the ROV survey methodology. We outline a spatial modeling approach that is capable of separating statewide and regional trends in the density of our focal species from effects attributable to the differences between MPAs and fished reference areas. We provide estimates of mean statewide (for four focal species with wide geographical ranges) and regional (north, central, and south California) trends in abundance and MPA effects across 11–17-year timeframes. Finally, we also compare the estimated MPA effect for each species with theoretical expectations based on data from fisheries stock assessments to explore whether these species are responding to protection as expected given their life histories and fishing pressure, thus providing a measure of MPA effectiveness in mesophotic depths across large scales for MPA management.

MATERIALS AND METHODS

Data collection and conditioning for analysis

ROV data collection and post processing methods were developed and tested by the California Department of Fish and Wildlife (CDFW) and Marine Applied Research and Exploration (MARE) and implemented in monitoring from 2005 onwards. Rectangular 500 m wide ROV survey sites were identified using high resolution seafloor maps and were placed perpendicular to the prevailing depth contour. MPA sites were paired with nearby reference sites in areas outside MPAs (Figure 1) that were chosen to cover similar habitats and depth profiles as much as possible (see Figure S1). Within MPAs and reference sites, the position for the shallowest 500-meter transect line was randomly generated and then additional transect lines were evenly spaced across the depth contour at the site (Figure 1). The number of transects selected at each site was based on the total percentage of rocky habitat present and the amount of transect effort needed to cover at least 4 km of that rocky habitat. A full description of data collection and conditioning methods used can be found in Lauermann et al. (2017). Collected video imagery was subsequently analyzed to characterize substrate types present and to identify and count all demersal and epibenthic finfish and macro-invertebrate species.

[IMAGE OMITTED. SEE PDF]

The ROV time-series data used in the present work spanned surveys across 26 MPAs and their associated reference areas that began in 2005 up until 2021, with between 2 and 9 repeat surveys (Table 1). The original ROV data set contained 3 MPAs that have been surveyed only once up until 2021, and two MPAs with 2 surveys only 1 year apart. Due to the uncertainty introduced in estimating trends through time with either a single survey or two surveys close in time, it was decided that these MPAs would be excluded from the present analysis.

TABLE 1 Summary of the time series of 500-m ROV transects used in the study across each marine protected area group (MPA group) including survey years, number of transects, total number of transects, and total number of repeats in each MPA group.

MPA group	Transects by year (500 m)
2005	2006	2007	2008	2009	2010	2011	2012	2013	2014	2015	2016	2017	2018	2019	2020	2021	Total transects	Total repeats
Point St. George Reef Offshore SMCA										23	14					19	12	68	4
Reading Rock SMR										19	19					20	14	72	4
Sea Lion Gulch SMR										15	6					18	20	59	4
Ten Mile SMR										19	20					20	18	77	4
Point Arena SMR/SMCA							12				17					14	12	55	4
Bodega Bay SMR/SMCA							31				45				38	44		158	4
Southeast Farallon Islands SMR/SMCA							21				27				23	23		94	4
Montara SMR											16					19	12.5	47.5	3
Pillar Point SMCA											8					12	9	29	3
Ano Nuevo SMR											9				10		10	29	3
Portuguese Ledge SMCA												15			12		10	37	3
Point Lobos SMR			12	31	23							24			23		34	147	5
Point Sur SMCA				22				25				23			22		20.5	112.5	4
Big Creek SMR/SMCA												28				13		41	2
Piedras Blancas SMR/SMCA												8				15		23	2
Point Buchon SMR				24	18			40				15			14		16	127	6
Campus Point SMCA										19					18		16	53	3
Harris Point SMR	30	24	21	21	19					23	23				24	33		218	9
Carrington Point SMR	25	31	25	25	25					25	26				24	40		246	9
South Point SMR	37	31	26	26	26					24	25				26	31		252	9
Gull Island SMR	44	41	39	39	38					39	40				32	41		353	9
Anacapa Island SMR/SMCA	39	29	30	28	25					59	29				28	37		304	9
Farnsworth Offshore SMCA										25						27	18	70	3
Swami's SMCA										25						13	14	52	3
Point Conception SMR										17					16		13	46	3
South La Jolla SMR/SMCA										24						27	20	71	3
Total	175	156	191	242	195		64	65		484	357	147			327	476	282	2841	119

MPAs across the California MPA network offer different levels of protection, with specific guidelines around activities allowed. The MPAs surveyed and included in the analyses were designated as either State Marine Reserves (SMRs) with full no-take protection, or State Marine Conservation Areas (SMCAs) with no-take and limited-take designations (Table 1). Prior to analysis, based on the specific rules for each included SMCA, it was determined that all surveyed SMCAs were no-take for all the focal species and therefore acted in the same way as SMRs in terms of direct protection, which was presumed to outweigh indirect impacts via prey or predator species in limited-take areas. SMRs and SMCAs were thus treated in the same way in assessing MPA effects in the analyses.

An analysis was performed for each species in turn. For four species with wide geographical distributions (copper rockfish, gopher rockfish, lingcod, and vermilion rockfish), we analyze the entire data set as “statewide” spatial models. For these four species, and an additional six species (brown rockfish, California sheephead, canary rockfish, kelp greenling, quillback rockfish, and yelloweye rockfish) we analyze trends using “regional” models. The regions used are those defined as long-term monitoring regions by the CDFW (CDFW & OPC, 2018): north (from the Oregon border to San Francisco Bay), central (San Francisco Bay to Point Conception), and south (Point Conception to the US/Mexico order). The form of all analyses was geostatistical spatial models (Lindgren & Rue, 2015), which in this instance involves describing each species abundance through space and time. The geostatistical model accounts for dependence between observations due to close proximity. Samples collected along transect lines, such as our subsampling units, are likely lack independence, and therefore methods such as our geostatistical model are required to deal with pseudoreplication and spatial autocorrelation (Elith & Leathwick, 2009). The variation in abundance is related to environmental descriptors, management (MPA) histories, as well as random spatial variation. Full details of the ROV data, the model structure, and the descriptors used (covariates) are given presently.

ROV data: Subunits for analysis

Associations between fish and preferred habitat are likely to be on a smaller scale than 500 m transects. Therefore, longer ROV transects are often split into subunits based on length or area for analysis. Previously, subunits ranging from 5m² (Enrichetti et al., 2023; Grinyó et al., 2018) to 50m² (Karpov et al., 2010; Karpov et al., 2012) in area have been used as well as subunits of 50 m (Duffy et al., 2014) and 20 m length (Budrick et al., 2019). The patchiness of habitats varied considerably amongst sites, but habitats were often patchy at scales of 10 m or less (see Figure S2), and it was reasoned that a sampling unit that captured the smaller scales of variation was preferable. Therefore, a subunit length of 10 m was chosen for analysis as it captures potential fine-scale habitat associations. Subunits of <6 m length (~ 6% of total) sometimes occurred when dividing longer transects up, and areas >50m² (~0.3% of total) could occur where the ROV was at a high altitude above the seafloor. Both were removed prior to analysis to avoid biases from short subunits or very large fields of view. A total of 133,506 10 m subunits were used in the final analyses.

Spatial model covariates: Environmental variables, temporal trends, and MPA effects

Both depth and habitat are known to be important drivers of the distribution of fish species. For example, many rockfish species are known to inhabit particular depth ranges, while most are associated with rocky reef habitat (Love et al., 2002). Habitat was included from the visual scoring of habitat classes in the ROV video recording. The start and end points of habitat classes (rock, sand, cobble, boulder, mud, and gravel) were scored continuously along transects. This allowed the subsequent calculation of the proportion of broad habitat classes (hard substrate, mixed substrate, and soft substrate) in each sampling subunit to be included as a habitat covariate. As proportions of habitat (hard, mixed, and soft) are constrained to sum to 1, including proportion soft habitat was redundant and was hence not included in the spatial model. To accommodate depth preferences, both a linear and quadratic effect function of depth were included in the spatial model. The quadratic term allows for “hump shaped” (i.e., “∩ shaped”) responses, which are expected when a species' preferred depth range is largely encompassed in the data.

Capturing large-scale variation in abundance using latitude alone is problematic for the California study region due to the nonlinear nature of the coastline, particularly moving south around Point Conception into the Southern California Bight. Therefore, for the statewide spatial models, to quantify large-scale spatial variation, a “coastal distance” variable was included. This variable represents the distance along a smooth representation of the 50-m contour (Figure 1). The smoothing was performed using a spline in the ArcGIS 10.1 software. This variable was set to zero in the north and extended to 1571 km at the southern-most surveyed site of South La Jolla. The resulting coastal distance variable was used in the model as a quadratic term (linear and squared terms as covariates), in a similar fashion to the depth² term (see above). For the regional spatial models, region (north, central, and south) was included as a categorical variable, as well as interaction terms between survey year and region and years since implementation (YSI) (see below) and region. This spatial model specification allowed the estimation of separate regional temporal trends and MPA effects. For the regional models, the coastal distance term was removed from the model as the region term and the interactions act as terms that capture large-scale variation at a regional scale and would otherwise be confounded with the coastal distance term.

A “survey year” term was included in all models to quantify the trend for any decrease/increase in density across the state or region not directly attributable to the MPA effect.

To quantify the temporal effect of MPAs separately from overall temporal trends experienced by both MPAs and reference areas, the cumulative effect of years of protection since MPA implementation was included as a covariate in the model (sensu Vanhatalo et al., 2017). We label this covariate as “YSI” and define it so that it is zero when an MPA is established and increases with time. Since the MPAs were not all established in the same year, the YSI variable varies across MPAs in any given survey year. All samples outside of MPAs (i.e., in reference areas) have zero years of implementation throughout the time series, whereas samples within MPAs were attributed the number of years since the MPA was established. Maximum values of YSI are located in the Channel Islands—the oldest MPAs, established in 2003 and last surveyed in 2020. The YSI term was included in the model using a log (YSI + 1) transformation to allow flexibility in how YSI impacts mean abundance in terms of the shape of the trajectory through time. Thus, for surveys in an MPA in the first year of implementation this formulation will set the cumulative effect at zero as log(0 + 1) = 0, and the effect in reference areas will remain zero throughout the time series. Modeling the temporal trends and MPA effects in this way means that the additional cumulative effect of MPA implementation is quantified in addition to any temporal trend shared by both MPA and reference sites.

The YSI effect (${\beta }_{\mathrm{MPA}}$) quantifies a power relationship (see Figure 2 for archetypical responses), where the estimate for YSI describes the trajectory in abundance inside MPAs compared with reference areas. In Figure 2, the red line results from a negative estimate for YSI and therefore a negative MPA effect indicating that MPAs are not achieving management goals of population increase compared with reference areas; the black line results from an estimate for YSI between 0 and 1, indicating in a positive response to MPA implementation that tends toward an asymptote; whereas the green line results from an estimate for YSI that is greater than 1 where there is a strong exponential growth inside MPAs compared with reference areas such as may occur in early stages of population recovery following heavy prior fishing pressure.

[IMAGE OMITTED. SEE PDF]

Description of the spatial modeling approach

All spatial modeling was conducted using integrated nested Laplace approximation (INLA; see Rue et al., 2009), which is a spatial regression approach that allows the incorporation of spatial effects, thereby accounting for spatial autocorrelation. These models use the locations ($s_i$) of the response variable (i.e., fish counts) and a vector of covariate values $X_i$ (that is ith row of the covariate matrix $X$) at those locations to estimate the expected values (${\mu }_i=E\left(y_i\right)$) at each location: 1 $\log \left({\mu }_i\right)=X_i\beta +\omega \left(s_i\right).$

The expected values (${\mu }_i$) depend on the covariate values ($X_i$) which are multiplied by coefficients ($\beta$'s) which are to be estimated, plus a spatial random effect $\omega$ that is estimated at each location ($s_i$). The spatial random effect is modeled as a multivariate normal distribution: 2 $\omega \left(s\right)~\mathrm{MVN}\left(0,\sum \right),$ with mean zero and variance–covariance matrix ($\sum$) that is populated using a Matérn correlation function (Matérn, 1986): 3 $\mathrm{cov}\left(i,j\right)=\delta \times \text{Matérn}\left(d_{\mathit{ij}},\kappa \right),$ where the covariance between any two locations $i$ and $j$ in the data set depend on their distance apart ($d$) in kilometers, the range of the Matérn function ($\kappa$), and the spatial variance ($\delta$) (for details about the INLA specifications see Lindgren et al., 2011). The default INLA smoothness setting for the Matérn function of α = 2 was used.

Thus, the spatial model parameters that need to be estimated are the coefficients of the covariate effects ($\beta$), the variation of the spatial effect ($\delta$), the range of the spatial effect ($\kappa$), and any residual variation. The spatial variance ($\delta$) quantifies the magnitude of spatially dependent variation, while the spatial range ($\kappa$) quantifies the distance that the correlation occurs over.

For the estimation of spatial random effects, a triangular INLA “mesh” that took into account the coastline and islands was defined across the state (see Figure S3). Sensitivity of results to the specification of the mesh was assessed with a sensitivity analysis (see Supporting Information).

The spatial model specification requires the specification of an appropriate error distribution. For these count data, we chose the negative binomial likelihood as fish count data can often exhibit overdispersion (relatively high variance with respect to the mean value). The area of each subunit was included as a model offset to account for the influence of different areas for each subunit, resulting in the response variable being modeled as density (i.e., number of fish per unit area).

A log link function is assumed to link the covariate effects to the expectation of the observations (Equation (1)). This model can then be expressed as 4 $\log \left(E\left(y\right)\right)={\beta }_0+{\beta }_{\mathrm{MPA}}\times \log \left(\mathrm{YSI}+1\right)+{\beta }_1{x}_1+\dots {\beta }_z{x}_z+\omega \left(s_i\right),$ where y is an observation, ${\beta }_0$ represents the model intercept and ${\beta }_1{x}_1+\dots {\beta }_z{x}_z$ represent other covariate effects, and $\omega \left(s_i\right)$ is the spatial random effect. Then, ignoring the other covariate effects (i.e., setting all coefficients to zero except for ${\beta }_0$ and ${\beta }_{\mathrm{MPA}}$): 5 $E\left(y\right)=e^{{\beta }_0}{e}^{{\beta }_{\mathrm{MPA}}\times \mathit{log}\left(\mathrm{YSI}+1\right)}={\left(\mathrm{YSI}+1\right)}^{{\beta }_{\mathrm{MPA}}}\times e^{{\beta }_0}.$

The exponentiated coefficient $\exp \left({\beta }_0\right)$ determines the mean density at the start of the survey ignoring all other covariates and spatial random effects. This specification means that for $0<{\beta }_{\mathrm{MPA}}<1$, the MPA effect will have a functional form with a shape that matches theoretical expectations outlined in Kaplan et al. (2019).

Spatial model priors

Bayesian models require the specification of prior distributions for parameters to be estimated. In particular, for a spatial model, the specification of priors for the spatial random effects (the spatial variance $\delta$, and the spatial range $\kappa$) are important considerations. We used “penalized complexity” priors (Fuglstad et al., 2019), and specified a prior probability of 0.1 that the spatial range was less than 2 km, and a prior probability of 0.1 that the spatial standard deviation was greater than 1. A sensitivity analysis was conducted by varying the priors an order of magnitude in either direction (i.e., spatial ranges of 0.2 and 20 km, and spatial standard deviations of 0.1 and 10) and exploring the impact on posterior estimates (see Figure S4), where it was determined that both these priors and the mesh used had minimal impact on inferences made.

The intercept determines the scale of the model, and its prior is specified for the precision (equal to the inverse of the variance), which was set at 0.0625. The priors for the precision of the remaining model coefficients were all set at 0.25. Environmental covariates of depth, proportion hard, proportion mixed, and coastal distance were scaled prior to modeling by centering with respect to means and scaling by the standard deviations.

Spatial model output summaries of temporal trends and the MPA effect

For both statewide and regional spatial models, temporal trends in abundance and changes in abundance within MPAs compared with reference areas were calculated and then visualized by taking 5000 joint posterior sample draws from fitted models. For each estimate, 95% credible intervals were calculated across the 5000 joint posterior draws. To avoid regional estimates with large uncertainties, regional temporal trends and MPA effects were only calculated where there was a minimum of 100 individuals of a species observed across the time series within a region. Thus, no regional assessment was made for brown rockfish in the south, California sheephead in the central or north, canary rockfish in the south, kelp greenling in the south, quillback rockfish in the central or south, or yelloweye rockfish in the south. When calculating temporal trends and MPA effects, mean environmental covariate values within the survey data were used and spatial effects were ignored (i.e., set to zero) as interest is in the specific fixed-effect response. Calculating the MPA effect and temporal trends in this way allowed either statewide or regional effects to be summarized as an average across the state/region, with estimates held at mean levels of the covariates.

Temporal trends were calculated as ratio changes in abundance through time. That is, the change in each year was calculated by using the estimated density in the relevant year in the state/region in the numerator and the estimated density in the initial survey year in the state/region in the denominator. Thus, for statewide models, ratios in temporal trends of abundance were calculated between 2005 and 2021; for the south region, they were calculated from 2005 to 2021; for the central region from 2007 to 2021; and for the north region from 2011 to 2021. For calculating the temporal trends for the statewide models, joint posterior draws of the intercept and survey year terms were used. For the regional models, joint posterior draws of the intercept (representing the abundance at the start of surveys in the central region), the regional fixed effects (representing starting differences to the central region for other regions), and the regional survey year estimates were used for estimating each regional trend separately.

For the MPA effect, we consider the ratio formed by the mean density after the YSI in the numerator versus mean density at the start of MPA implementation using Equation (5). MPA effect ratios were calculated over 17 years of protection for statewide models, 17 years of protection for the south region, 14 years of protection for the central region, and 11 years of protection in the north for all regional models. This construction allows for direct comparison with theoretical expectations developed under an alternative predictive population dynamics model (see “Comparison between spatial model estimates and theoretical MPA expectations from population dynamics models” section below).

Comparison between spatial model estimates and theoretical MPA expectations from population dynamics models

The estimated MPA effect quantified by the YSI coefficient was compared with theoretical MPA responses for each species outlined in Kaplan et al. (2019). Theoretical expectations were based on simulations from age-structured population models taking into account life history parameters of each species, and recruitment variability and fishing mortality estimates from stock assessments. Demographically open population (i.e., high larval immigration) models were used. Full details and parameters used can be found in Kaplan et al. (2019). Note that Kaplan et al. (2019) did not provide predictions for canary rockfish, quillback rockfish, or yelloweye rockfish, and so for these species no comparison could be made. Theoretical expectations and the 95% confidence intervals from the simulations were overlaid with plots of the estimated MPA effects and 95% credible intervals from our models to allow direct comparisons.

The ROV data and all code used for modeling and plotting are available at: .

RESULTS

Temporal trends in the density of focal species

For the statewide spatial models, the coefficients for survey year were positive for all species (Figure 3 and Supporting Information), indicating an overall increase in density for copper rockfish, gopher rockfish, lingcod and vermilion rockfish statewide between 2005 and 2021. Gopher rockfish showed the largest increases with an estimated approximate 17-fold increase between 2005 and 2021, followed by copper rockfish with an approximate 6-fold increase, vermilion rockfish with a 3-fold increase, and lingcod with a 2-fold increase (Figure 3).

[IMAGE OMITTED. SEE PDF]

For the regional models, the coefficients for survey year were positive for 18 out of 22 species-region combinations assessed, no significant effect for four species-region combinations, and only one negative estimate (lingcod in the north region) (Figure 4 and Supporting Information). Non-significant estimates for temporal trends (i.e., 95% credible intervals overlapping zero) were found for canary rockfish in the central region, and kelp greenling in all regions. Positive mean regional estimates of abundance increase ranged from modest levels of approximately 1.5-fold for lingcod in the central region 2007–2021, to increases of greater than 10-fold for brown rockfish and copper rockfish in the central region 2007–2021, and gopher rockfish in the south region 2005–2021 (Figure 4). In some instances, there was considerable variation in temporal trends for the same species between the different regions. For example, gopher rockfish showed higher trajectories in the south and central regions than the north; lingcod showed higher trajectory in the south region, a low positive trajectory in the central region, and a negative trajectory in the north region; and vermilion rockfish showed a higher trajectory in the central and north regions compared with the south (Figure 4). Credible intervals were also notably wider in the regional models compared with the statewide models (Figures 3 and 4).

[IMAGE OMITTED. SEE PDF]

MPA effects for focal species: Statewide and regional estimates

For the statewide spatial models, positive coefficients for YSI were found for all four species with all 95% credible intervals not containing zero (Figure 5 and Supporting Information), indicating the positive influence of MPAs on the density of all modeled species. The strongest responses were for copper rockfish and gopher rockfish with an estimated approximate 2.9-fold and 2.3-fold increase in density inside MPAs over the 17 years of MPA protection (Figure 5). Lingcod and vermilion rockfish both displayed smaller responses of ~1.5 and 1.3-fold mean increases over the 17 years of MPA protection, respectively.

[IMAGE OMITTED. SEE PDF]

For the regional models, positive coefficients for YSI were found for 18 out of 22 of the species/region combinations modeled, with four species/regions combinations displaying non-significant YSI coefficients, and no species/region combination displaying negative coefficients (Figure 6 and Supporting Information). Non-significant MPA effects (i.e., 95% credible intervals containing zero) were found for gopher rockfish and vermilion rockfish in the south region; brown rockfish, yelloweye rockfish in the central region; and yelloweye rockfish in the north region. Estimates of regional abundance changes ranged from relatively small (e.g., an ~1.5-fold increase for lingcod in the central region) to greater than 10-fold increases (e.g., brown and copper rockfish in the central region and gopher rockfish in the south region). Credible intervals were again wider for regional estimates of MPA effects compared with statewide estimates (Figures 5 and 6).

[IMAGE OMITTED. SEE PDF]

Comparisons of the estimated MPA effect with theoretical expectations

For the statewide models, all species showed a positive increasing MPA effect, with estimates for 0 < YSI < 1 (Figure 5 and Supporting Information) and therefore followed a trend similar to that of the black line in Figure 2. Mean MPA responses estimated from the ROV data for both copper and gopher rockfish exceeded the theoretical expectations over the 17 years of protection (Figure 5). Mean responses estimated from the ROV data were lower than expected for both lingcod and vermilion rockfish, although there was still considerable overlap in the credible/confidence intervals between the spatial models and population dynamics models of the effects for both species.

For the regional models, 11 out of 17 species-region comparisons showed higher increases in abundance than expected. These were California sheephead, copper rockfish, and lingcod in the south; brown rockfish, copper rockfish, gopher rockfish, kelp greenling, lingcod, and vermilion rockfish in the central region; and brown rockfish, copper rockfish, and gopher rockfish in the north region (Figure 6). Species that showed lower than expected responses were lingcod in the south region, and kelp greenling, lingcod, and vermilion rockfish in the north region. In some cases, regional estimates of MPA effects were quite different between regions. For example, increases inside MPAs were relatively high and above expected in the central region for lingcod (~4-fold), but below expected in the south and north regions (both ~1.8-fold). For vermilion rockfish, estimates were ~ 2.2-fold and ~ 1.6-fold for the central and north regions, respectively, but ~0.8-fold (i.e., slightly decreasing, but not statistically significant—see Supporting Information) in the south region.

Additional covariate and spatial effects

Positive effects were found for all species for the proportion of hard bottom habitat and proportion of mixed habitat from the visual ROV data (Supporting Information). In the statewide spatial models, positive depth effects were found for copper rockfish, lingcod, and vermilion rockfish indicating an increase in density with depth, and a negative effect for depth was found for gopher rockfish indicating a preference for shallower depths across those surveyed (Supporting Information). Expected habitat preferences for rocky reef were also found for the focal species in the regional models, with brown rockfish, canary rockfish, quillback rockfish, and yelloweye rockfish showing preference for deeper surveyed depths, while California sheephead and kelp greenling showed preferences for shallower depths (Supporting Information).

Spatial dependence (i.e., a correlation in abundance in space) was found to exist across scales of 1.6 km (vermilion rockfish) to 5.3 km (lingcod) across the four species in the statewide spatial models, and 1.7 km (vermilion rockfish) to 507 km (yelloweye rockfish) in the regional spatial models (range estimates, Supporting Information). The spatial standard deviation was between 1.7 (copper rockfish) to 2.5 (gopher rockfish) in the statewide spatial models, and 1.1 (kelp greenling) to 9.4 (brown rockfish) in the regional spatial models. This level of spatial variation is not insignificant and occurs over relatively small scales, with the exception of yelloweye rockfish, which is generally low abundance over the coastline, indicating that abundance of the majority of the focal species is patchy along the Californian coastline.

DISCUSSION

Here, we present results of the analysis of a 17-year (2005–2021) ROV monitoring program across California's MPA network, highlighting increases in statewide and regional abundance for 10 focal fish species, as well as the effectiveness of the California MPA network in promoting additional increases in abundance within MPAs. Our findings of MPA effectiveness across mesophotic depths in California's MPA network support previous findings of MPA effectiveness in shallower depths (Caselle et al., 2015; Hamilton et al., 2010), while our findings of increases outside MPAs are supported by many of the recent stock assessments for the species we assessed (NOAA Fisheries, 2024). At a statewide-level, we observed overall increases in the abundance of all four assessed species (copper rockfish, gopher rockfish, lingcod, and vermilion rockfish) between 2005 and 2021, with abundance levels increasing between 2 and 17-fold on average. Regional models (based on south, central, and north long-term monitoring regions) also detected strong recovery in abundance within each region, with 18 out of 22 region-species combinations demonstrating increases in abundance, three showing non-significant changes, and only one (lingcod in the north region) showing a decline in abundance. We show that MPAs contributed additional abundance increases compared with fished areas statewide for all 4 species assessed, and for 10 species assessed on a regional basis, 18 out of 22 species-regions assessed displayed positive MPA effects, 4 showed no detectable effects, and none showed negative MPA effects. Our findings with respect to MPA effects when compared with theoretical expectations from population dynamics models showed that 2 out of 4 species across the statewide MPA network demonstrated recovery rates that exceed what was anticipated within MPAs, and 11 out of 17 exceeded expectations within MPA monitoring regions. These results have important implications for management, as they provide evidence of the California MPA network's ability to meet goals for the restoration of natural abundance of previously fished species under the Marine Life Protection Act (1999) and indicate that there has been significant statewide and regional recovery in our focal species density over time. This information is particularly useful for the adaptive management of the MPA network, and more broadly for informing fisheries management of these species outside of the MPAs. Our findings also highlight the usefulness of imagery-based surveys for assessing MPA network effectiveness across mesophotic depths and tracking the abundance of species that are of interest for fisheries management purposes over time.

We found strong statewide and regional trajectories of increased density for all four species at a statewide-scale between 2005 and 2021, and for 18 out of 22 species-region combination at a regional scale over 10–17 years of monitoring. These findings are in general agreeance with recent stock assessments for many of these species, with predicted increased spawning stock biomass across California for copper rockfish (Wetzel et al., 2021a, 2021b) and vermilion rockfish (Dick et al., 2021; Monk et al., 2021) since the early 2000s, and lingcod since approximately 2010 in the south (Johnson et al., 2021), and since the 1990s in the north (Taylor et al., 2021); and in wider stock assessments across the US Pacific coast for brown rockfish between 2000 and 2013, canary rockfish between 2000 and 2023, and yelloweye rockfish between 2001 and 2017 (NOAA Fisheries, 2024). Two notable exceptions are gopher rockfish and quillback rockfish, which are both noted to have declines in recent California stock assessments. The 2019 stock assessment for gopher rockfish showed a predicted decline in spawning stock biomass since the early 2000s (Monk & He, 2019) in contrast to our predicted large increase in abundance. However, another monitoring program based on catch data also reports large increases in abundance of gopher rockfish since 2013 (Hamilton et al., 2021). Quillback rockfish also show a decline in estimated stock abundance across the California coast between 2009 and 2021 (NOAA Fisheries, 2024) in contrast to our predicted increase. Stock assessments are generally reliant on catch data and therefore may not capture recent pulses of new recruits that are not contributing to spawning stock biomass estimates but that can be quantified in ROV video footage. While the drivers of these outside trajectories were not the focus of this study, successful recruitment years preceding and over the survey period for these species as well fisheries management measures such as quota limits and depth restrictions (i.e., rockfish conservation areas) are likely to contribute positively to abundances. The MPA network is also expected to benefit outside areas through “recruitment subsidies” (e.g., Le Port et al., 2017) via increased larval export, and potential “spillover” (e.g., Di Lorenzo et al., 2020) where fish migrate outside MPA boundaries. We do not directly quantify spillover effects; however, our spatial modeling approach quantified spatial dependence for different species on scales of 3–8 km, which often includes nearby reference areas, and thus may be capturing some localized spillover. More detailed studies such as those that utilize knowledge about home ranges of species, MPA-specific habitat fragmentation and outside fishing pressure (Pinillos & Riera, 2022), or detailed genetic and/or biophysical modeling studies (e.g. Christie et al., 2010) may provide further evidence of spillover effects or larval export from MPAs to outside areas.

One of the expected responses to MPA establishment is the increased density of previously fished species through time (Claudet et al., 2008; Hopf et al., 2016). We found that all four species assessed at a statewide-level, and 18 out of 22 species-region combinations assessed displayed positive trends in density inside MPAs over 11–17 years of protection. A meta-analysis of 80 MPAs found that density increases in previously fished species were often detectable in 3 years or less (Halpern & Warner, 2002). However, many of the species in the present study are long-lived (Love et al., 2002) with sporadic boom-bust recruitment dynamics, which may result in delayed detectability of MPA effects (Hopf et al., 2021). The positive MPA effects found for the vast majority of species reported here highlights that the California MPA network has been in place for a sufficient length of time for expected increases in density to be detectable for these species.

The timely reporting of protection effects is critical to management of MPAs and is reliant on well-designed monitoring programs that are capable of capturing trends through space and time. The ROV surveys used in this study spanned 17 years and over 1500 km of coastline covering mesophotic depths (~20–150 m), depths which encompass the majority of reef area protected by California's MPAs (Starr et al., 2021). Previous work has shown that detecting MPA effects at the scale of individual MPAs in the California study system is likely to require high levels of within-site sampling and long periods (>15–20 years) of data collection (Perkins et al., 2020; Starr et al., 2015). For the analyses presented here, data across 26 MPAs and associated reference sites was used, with the spatial modeling approach able to capture average statewide and/or regional trends and MPA effects while accounting for the differences between MPAs that was not captured by the environmental covariates. The scope of this data and the modeling approach used allowed for estimates of temporal trends and MPA effects to be made for each species at network-wide and management region scales. This is especially important for networks of MPAs, such as California's MPA network, where the imperative is to report and manage the MPAs on a network or bioregional scale (Hall-Arber et al., 2021).

The comparison of observed rates of recovery due to MPA establishment with theoretical expectations is an important component of an adaptive management framework (Kaplan et al., 2019; Nickols et al., 2019). Our results show that on a statewide scale, expectations are being exceeded for copper and gopher rockfish, and on the lower range of expectations for lingcod and vermilion rockfish. For our regional assessments, 11 out of 18 species-region combinations exceeded expectations including California sheephead and copper rockfish in the south region; brown rockfish, copper rockfish, gopher rockfish, kelp greenling, lingcod, and vermilion rockfish in the central region; and brown rockfish, copper rockfish, and gopher rockfish in the north region. As well as life history characteristics for each species, theoretical expectations used estimates of fishing mortality and recruitment variability based on stock assessments and averaged over the years prior to MPA implementation (Kaplan et al., 2019), both of which are likely to be spatially variable. The timing and magnitude of recruitment pulses are likely to have a large impact on the population dynamics and detectability of MPA effects (Hopf et al., 2021). This suggests that further exploration of recent recruitment dynamics for these species as well as the effects of spatially variable fishing mortality may be beneficial in understanding the results presented here. Furthermore, we note that the theoretical expectations presented are based on life history characteristics that may in some cases be poorly characterized, despite the spatial averaging used in stock assessments, which adds additional uncertainty to these theoretical estimates. However, these are the current best estimates of expected trajectories and thus are the best available information for present comparisons.

Our spatial modeling approach found moderate to high spatial correlation in abundance at local scales (several km) as well as differing temporal trends and MPA effects for the same species between regions. For example, vermilion rockfish were found to have a lower-than-expected MPA effect in the statewide spatial models that seems to be primarily driven by trends in the south region, where both temporal trends and MPA effects were lower than the other regions. This indicates quite different dynamics in this region for vermilion rockfish that warrants further exploration. Estimates of the posterior spatial range and spatial standard deviation from our models indicate a high level of variability in density at the scale of individual MPAs and their associated reference areas, indicating localized (km scale) processes are important drivers in the density of the species modeled at scales smaller than regions. Modeling residual spatial autocorrelation that is not accounted for by model covariates has been shown to be important, as coefficient estimates (such as YSI) can otherwise be erroneous (Dormann, 2007). Spatial autocorrelation across smaller scales is likely to be particularly important for data collected along transects such as the ROV data analyzed here. High spatial variability in MPA responses has been noted in SCUBA diver monitoring data in the Channel Islands (Caselle et al., 2015) where incorporating biogeographical differences helped explain some of these differences (Hamilton et al., 2010). While full exploration of potential drivers of spatial variability in responses is outside of the scope of the present work, we note that our spatial modeling approach captures unaccounted for spatial variability in the spatial random effects and thus offers a preferred approach for analyzing such data sets where high spatial variability is likely to be present.

Here, we demonstrate positive statewide and regional increases in abundance of 10 focal species as well as additional positive MPA effects across mesophotic depths in California's MPA network using video-based ROV surveys. Our results support many recent stock assessments that report increasing trends in abundance, as well as previous reports of MPA effectiveness across shallower depths. Quantifying the effects of management interventions, such as MPAs, in the marine environment is challenging due to the technical difficulties in collecting data, particularly in mesophotic depths. Our results demonstrate that technological solutions, such as ROVs, allow collection of data across spatial scales, depths, and timeframes capable of providing data that can quantify MPA effectiveness for MPA managers as well as inform wider trends in abundance for fisheries management.

AUTHOR CONTRIBUTIONS

A.L. and M.P. collected the data and contributed to writing and editing the manuscript. G.R.H. and S.D.F. assisted with the statistical analysis and writing and editing the manuscript. N.R.P. conducted the analyses and wrote the initial draft.

ACKNOWLEDGMENTS

The authors would like to acknowledge the California Ocean Protection Council that funded some of the ROV survey work and subsequent annotation and analysis. We would also like to thank CDFW for establishing the ROV monitoring program and for their ongoing support and feedback. We thank Rosie Byers and Nissa Kreidler for their help with the GIS and mapping component of the project. We would also like to thank all crew members involved in the deployment of the ROV over the years. Open access publishing facilitated by University of Tasmania, as part of the Wiley - University of Tasmania agreement via the Council of Australian University Librarians.

DATA AVAILABILITY STATEMENT

The data used in analysis and code are available at .

References

Agardy, T. (2000). Information needs for marine protected areas: Scientific and societal. Bulletin of Marine Science, 66(3), 875–888.

Budrick, J., Ryley, L., & Prall, M. (2019). Methods for using remotely operated vehicle survey data in assessment of nearshore groundfish stocks along the California coast. Retrieved from https://www.pcouncil.org/documents/ on March 19th 2020.

Caselle, J. E., Rassweiler, A., Hamilton, S. L., & Warner, R. R. (2015). Recovery trajectories of kelp forest animals are rapid yet spatially variable across a network of temperate marine protected areas. Scientific Reports, 5, [eLocator: 14102]. [DOI: https://dx.doi.org/10.1038/srep14102]

Wetzel, C. R., Langseth, B. J., Cope, J. M., & Budrick, J. E. (2021a). The status of copper rockfish (Sebastes caurinus) in U.S. waters off the coast of California north of Point Conception in 2021 using catch and length data. Portland, Oregon.

Wetzel, C. R., Langseth, B. J., Cope, J. M., & Budrick, J. E. (2021b). The status of copper rockfish (Sebastes caurinus) in U.S. waters off the coast of California south of Point Conception in 2021 using catch and length data. Portland, Oregon.