We replicated the sampling methods of Rudershausen et al. (2010) and collected data inside the SWMPA along the continental shelf break in 47–150 m deep. Besides the shipwreck near the offshore edge of the area, the only known habitat for reef fish concentrations is along the inshore edge of the SWMPA (Rudershausen et al., 2010). Our sampling was therefore restricted to a 4.63 × 27.78 km rectangle inside the SWMPA as well as an adjacent area of equal dimensions that has never been closed to bottom fishing (hereafter: “control” or “control area”; Figure 1) that was also sampled in Rudershausen et al. (2010). Data collection trips from 2018 to 2020 were attended by at least one author of Rudershausen et al. (2010) and/or the collaborating boat captain from that study.

The null hypothesis in this study was that the SWMPA had no effect on overall biomass, abundance, individual size, and individual age of reef fishes. Within each area, we collected hydroacoustic backscatter and biological data to test this hypothesis. Hydroacoustic backscatter data were collected along a zig-zag track that connected alternating points along the inshore and offshore long edge of each rectangular area. Points along the edges of the areas through which we transected were spaced 4.63 km apart and specific paths were varied between trips so that we sampled new habitat each day. This followed the sampling protocol for the majority of the Rudershausen et al. (2010) data collection, though in 2007 the authors of that study conducted some transects in the SWMPA both parallel and perpendicular to the long edge of the area. For a given day of sampling, data were collected in either the SWMPA or the control area.

The modern study used an echosounder consisting of a Simrad ES60 transceiver outfitted with a single-beam transducer that operated at 38 kHz. We elected to use this gear because it was the exact unit used in Rudershausen et al. (2010). Acoustic data were collected from the R/V Cape Fear (20-m diesel-powered former fishing vessel) in 2018–2019 and from the R/V Regulator (8-m center-console vessel) in 2020. On both vessels, the transducer was mounted to the base of an aluminum pole and deployed off the port side so that the face of the transducer was approximately 0.5 m below the water surface. Boat speed ranged from three to five knots during acoustic data collection. The transceiver was connected to a laptop computer which was also connected to an on-board GPS plotter for positioning. During sampling, we noted the presence of any schools of on- or near-bottom (within 10 m of the seafloor) biomass (“acoustic events”) and recorded their GPS locations. Raw ES60 files were saved for subsequent processing. These methodologies are consistent with those of Rudershausen et al. (2010).

A subset of acoustic events were biologically sampled with hook-and-line. We determined ad hoc which events would be sampled with which gear based on a combination of the observed size of each acoustic event (to sample a range of fish group sizes) and their proximity to our vessel’s location (for logistical reasons). The choice of which acoustic events to sample was consistent between the SWMPA and control areas. Hook-and-line sampling gear consisted of conventional rods with electric reels and 59 kg braided line. Terminal tackle was a high-low bottom rig made of 68 kg monofilament line and two 8/0 J-style hooks and lead weight ranging from 0.68 to 1.36 kg. Hooks were baited with cut squid (Illex or Loligo sp.). Hook-and-line sampling usually took place within 24 h of acoustic data collection, though in one instance a hook-and-line trip took place 23 days after acoustic data collection.

When hook-and-line fishing, captains were permitted to keep vessels in gear (i.e., hover over events) or out of gear (i.e., drift) depending on sea conditions. Gear was deployed within ∼50 m of an acoustic event. Hook-and-line sampling was terminated when the boat drifted farther than ∼50 m from the event, whereupon the captain repositioned the boat. Hook-and-line sampling continued at an event until at least four drops of baited rigs had been conducted; a drop consisted of one vertical round-trip of a two-hook rig that was fully baited when deployed. Catch per unit effort (CPUE) was measured as fish-per-drop (regardless of species) for each baited two-hook rig. Caught reef fish were identified, measured (fork and total length; mm), and released. In 2019–2020, red porgy (Pagrus pagrus) were retained for aging.

Acoustic data were processed using Echoview software (v. 10.0.275, Echoview Party Ltd.). We processed data that were collected during the present study (2018–2020) as well as raw acoustic data from the Rudershausen et al. (2010) study (collected in 2007–2009). For each data file, seafloor definitions were edited by one or two authors (BR and/or EE) and an impulse noise removal filter was applied to all files to reduce systematic non-biomass backscatter. Data were limited to backscatter stronger than -60 dB. For one analysis, we imposed a 100 m linear grid on each data file. Grid cells were bounded vertically by the seafloor and a line 10 m above it. This analysis yielded a total backscatter value for each 100 m × 10 m cell as the Nautical Area Scattering Coefficient (NASC; also called S_A; units m² nmi^–2). The Echoview software calculated ensonified water volume accounting for beam spreading and scaled NASC values depending on depth. Additionally, we expected that habitat quality would be correlated with biomass; the continental shelf break inside the SWMPA and control area consists of stretches of low-slope unconsolidated sediment punctuated by high-relief drop-offs and ledges (Quattrini and Ross, 2006; Rudershausen et al., 2010). High vertical relief habitat is known to aggregate greater biomass of demersal fishes (Randall and Minns, 2000; Claisse et al., 2014), therefore we wished to include this metric in our statistical model. We generated an estimate of seafloor slope for each 100 m grid cell in our survey. Echoview exports included mean depth for each 100 m cell, so we estimated slope by generating a raster image for the inshore edge of the SWMPA and the control area using all mean depth observations from all acoustic sampling trips. Point values for mean depth were used to interpolate full coverage of each area using inverse distance weighting. Slope in degrees for each cell was then calculated using a 3 × 3 cell moving window of mean depth. Because this methodology relies on depth information from adjacent locations, we could not use it to estimate slope for some grid cells around the edge of each area. For these locations, we applied the slope from the nearest (linear distance) grid cell in which it was estimable. These analyses were conducted in R (R Core Team, 2021) and Esri ArcPro (v 2.7).

For a separate analysis, individual acoustic events (aggregations of fish) were defined and their acoustic backscatter was measured (again in terms of NASC). Only acoustic events that were sampled biologically in our study or during Rudershausen et al. (2010) were measured for this second analysis. The definition of what was considered a fish aggregation, determined by a single author (EE), was consistent across data files (and thus between areas and time periods), and made no assumptions about species composition. In general, fish aggregations were delineated by examining the echogram for clear areas of higher density; when several areas of higher density were adjacent to one another, they were usually considered distinct aggregations (depending on their size and spatial distribution).

We fitted a series of lognormal hurdle models in a Bayesian framework to the 100 m grid cell hydroacoustic data using the R package “brms” (Bürkner, 2017). The first stage of the hurdle model estimated the probability of a given cell containing zero biomass and the second stage estimated the distribution of biomass given that it was non-zero. We fitted a range of candidate models to the data. The full model included Area (SWMPA or control) and Period (before or after) as main effects as well as their two-way interaction. In a BACI design such as our own, an Area:Period interaction tests for an MPA effect; if the directionality and/or magnitude of change through time differs between the two areas, this suggests that designation of the MPA (or “impact”) is the cause. In addition to these main effects, we included a fixed effect for seafloor slope (Slope; continuous), a fixed effect for season (Season; factor; winter, spring, summer, and fall), and a random effect for trip (Trip; intended to account for stochasticity in environmental conditions that we were unable to measure). We considered including seafloor depth as a variable in these models, however, ad hoc investigations revealed that depth and slope were highly correlated. This was unsurprising considering the most extreme drop-off of the continental shelf tends to occur around the same depth regardless of latitude. The correlation between these two variables was non-linear (Supplementary Figure 1), and a generalized additive model (GAM; R package “mgcv”; Wood, 2011) revealed that the relationship was highly significant (Supplementary Figure 2; edf = 8.5; p < 0.001). We therefore excluded depth as a predictor variable in all models.

The full model was specified as:

Where hu_β₀ through hu_β₅ and β₀ through β₅ were given a prior of Normal(0,10), hu_α_j and α_j were given a prior of Normal(0, hu_σ_Trip) and Normal(0, σ_Trip) respectively, and hu_σ_Trip and σ_Trip were given a prior of HalfCauchy(10). These priors were generally uninformative but may have weakly constrained the parameter estimates to exclude outlying observations. In this model, Y_i is the amount of biomass (NASC) per 100-m grid cell in the sonar survey, β terms describe coefficients for each variable, α_j is a random intercept for trip j, and σ_j is a random standard deviation for the intercept term for trip j. Other candidate models included different combinations of these variables, starting with the most basic model (intercept only) and working up to the full model. Candidate models were examined only if they were deemed biologically plausible and no three-way interaction terms were included in any model (Table 1). All models used the same effects for both model stages, as we assume that a covariate that is associated with high encounter probability will also likely be associated with high biomass (Thorson, 2018). Each model was fitted assuming a lognormal distribution and using four sampling chains, each with 10,000 iterations, a burn-in period of 1,000 iterations, and thinning of every third sample. We examined model output, including R-hat and effective sample size, to ensure convergence. We used Widely Applicable Information Criterion (WAIC) to compare our candidate models and select the best model; WAIC is useful for comparing complex Bayesian models and often outperforms other metrics such as AIC or DIC (Watanabe and Opper, 2010; Gelman et al., 2014). Meaningfulness of partial regression coefficients was evaluated by examining 95% credible intervals and checking whether they contained zero.

We also fitted hurdle models to the data consisting of Catch Per Unit Effort (CPUE) at sites that were biologically sampled. In this procedure, the first stage of the hurdle model estimated the probability of a given site producing zero CPUE (no fish caught) and the second stage estimated the distribution of CPUE given that it was non-zero. We again fitted a range of candidate models to the data, beginning with an intercept only model and working up to various combinations of the main effects (Area and Period) and the Area:Period interaction along with a fixed effect of biomass (NASC) calculated at each individual site (Table 2). Finally, we included a random effect of Trip in some candidate models, intended to account for among-day variability in sea conditions and fish catchability.

The full model was specified as:

All models used the same effects for both model stages, and all models were fitted assuming a lognormal distribution using four sampling chains, each with 10,000 iterations, a burn-in period of 1,000 iterations, and thinning of every third sample. Prior distributions, model evaluation, and parameter evaluation occurred as above. We replicated this modeling procedure for a subset of CPUE data consisting of only species that have formal stock assessments in this region, which is an indication of their economic importance and a reasonable proxy for how heavily they have been targeted in recent decades. For this analysis, we included CPUE data for only blueline tilefish (Caulolatilus microps), gag (Mycteroperca microlepis), gray triggerfish (Balistes capriscus), greater amberjack (Seriola dumerili), red grouper (E. morio), red porgy (Pagrus pagrus), scamp (M. phenax), snowy grouper (H. niveatus), and vermilion snapper. Finally, we performed a qualitative investigation of overall CPUE trends in the four crosses of Area and Period for CPUE of all reef fish and CPUE of assessed species only. This analysis was conducted with per-trip CPUE.

We examined reef fish community composition in both periods and areas using non-metric multidimensional scaling (NMDS) and analysis of similarities (ANOSIM). NMDS and ANOSIM have been used to characterize fish communities in impacted areas as compared to control sites (Shepherd et al., 1992; Bacheler et al., 2016). We treated each area/period cross as an independent community of reef fish. Given that gear was standardized between time periods, we assumed that changes in catch reflected changes in the community rather than an effect of sampling. Our sampling unit was each trip during which reef fish were caught; we could have divided the data to the site level, but many sites produced zero catch which would require their exclusion from community analyses. We binned reef fish into six groups: groupers (Serranidae), jacks (Carangidae), porgies (Sparidae), snappers (Lutjanidae), tilefish (Malacanthidae), and triggerfish (Balistidae) (Table 3). Observations of species that did not fit into these groups were rare (<5% of all catch) and were therefore excluded from the community analyses. Catch data were standardized by effort (number of drops) for a given trip, yielding a CPUE value for each trip and species group; relative rather than absolute values are required for NMDS when sampling effort is variable (Clarke, 1993). CPUE data were fourth-root transformed to down-weight the dominance of highly abundant species groups (Clarke and Warwick, 2001). Relationships between crosses of area and period were visualized by plotting NMDS results based on zero-adjusted Bray-Curtis similarity; this adjustment controls for erratic behavior that sometimes occurs for zero-rich data (Clarke et al., 2006). We conducted a one-way ANOSIM using the same zero-adjusted Bray-Curtis similarities as in the NMDS analysis. The ANOSIM was used to calculate an R statistic and differences between the four communities were evaluated for significance at α = 0.05. Multivariate analyses were performed in the R package “vegan” (Oksanen et al., 2020).

We conducted single-species analyses with red porgy for several reasons. First, they were the most frequently caught species in both areas and time periods; other species lacked appropriate sample sizes in one or more areas or periods (see “Results”). Moreover, red porgy have fairly narrow home ranges relative to the size of the SWMPA; short term ranges are generally ∼0.5 km² (Afonso et al., 2009) while the entire SWMPA is >500 km². Red porgy reach reproductive maturity quickly, with most spawning by age 2 or 3 (Manooch and Huntsman, 1977; Wyanski et al., 2019). Given the size and age of the SWMPA relative to these biological characteristics of red porgy, this species represented an ideal candidate for testing a monospecific MPA effect. Moreover, for similar reasons, Pickens et al. (2021) conducted single-species analyses for red porgy caught in three of the other SEUSA MPAs; by singling out this species, we added to their analysis of how SEUSA MPAs are affecting biological metrics of red porgy.

For red porgy collected in 2019–2020, fish were weighed (whole weight; g) and sagittal otoliths were removed. Otoliths were processed and analyzed at the NOAA Fisheries Laboratory in Beaufort, North Carolina, using the methods described by Burton et al. (2017). Calendar ages were assigned to each specimen based on the number of opaque annuli and the month of capture (J. Potts, NOAA Fisheries, personal communication). In addition to red porgy collected by our survey, we obtained age information for red porgy caught by the Southeast Reef Fish Survey (SERFS) in the SWMPA and the surrounding region of Onslow Bay, NC, including the control area, from 2015 to 2019 (T. Smart, SC DNR, personal communication). SERFS data were only available through 2019, so red porgy encountered in 2020 for both inside control and outside SWMPA were all from our collections in the control area. Sample sizes in the SWMPA prior to 2015 were insufficient for inclusion. Similarly, no red porgy ages were available from Rudershausen et al. (2010), which precluded a BACI framework for this analysis. Red porgy collected in the present study were collected with hook-and-line and chevron traps (B. Runde, unpublished data); red porgy from SERFS were collected with chevron traps and short bottom longlines. We assumed that the selectivity of these gears is equivalent because a recent stock assessment of this species used logistic selectivity curves for hook-based gear as well as traps, and both curves used age 6 as the youngest age of full recruitment to the gear (SEDAR, 2020). We examined age frequencies of red porgy by area (SWMPA or control) and year (2015–2020). We calculated mean age by area and year and compared yearly means using a Wilcoxon rank-sum test. Significance was assessed at α = 0.05 with a Bonferroni correction applied such that α = 0.008. In addition to mean age, we compared yearly age diversity between the SWMPA and control sites using the Shannon diversity index. Age diversity is an important metric of fish stock health, with higher diversity implying greater health (Marteinsdottir and Thorarinsson, 1998). Diversity values were calculated using the R package “vegan” (Oksanen et al., 2020).

We compared size (TL, mm) distributions of red porgy caught by the authors of Rudershausen et al. (2010) with those caught in 2018–2020 in both areas. To augment our dataset in the latter period, we obtained length information for red porgy caught by SERFS in these two areas (T. Smart, personal communication). As in the age analysis, we assume selectivity of gears used in these surveys was equivalent for red porgy. For the size analysis, we did not use red porgy data from elsewhere in the region; analyses were restricted to individuals that were caught in either the SWMPA or the control area. We did not expand our dataset to include fish from outside the control area because data were not available from SERFS in the “before” period. We used Bayesian linear model to compare the four size distributions (TL, centered) and tested for the effect of Area, Period, and their interaction. Meaningfulness was evaluated based on whether or not the credible intervals contained zero.

Across both areas and periods, sonar data collection resulted in 10,951 grid cells of 100 m length. The distribution of non-zero density (NASC) observations in these cells was approximately lognormal (Figure 2). The best model (lowest WAIC) was Period + Slope + Season + Trip (Table 1). There were several additional models that were competitive in terms of WAIC; each of these models contained Slope and Trip. Results from the best model indicated meaningful effects of Slope in both stages of the model; cells with higher slope were less likely to contain zero biomass, and given non-zero biomass, cells with higher slope had higher biomass (Table 4 and Supplementary Figure 3). In this model, the only other variable whose credible interval for its partial regression coefficient that did not include zero was Period, for which the before period had less biomass than the after period (Supplementary Figure 4). We also investigated results from the full model given that it was nearly equivalent to the best model in terms of WAIC (ΔWAIC = 0.5, Table 1). The only variable whose partial regression coefficient did not include zero in either stage of this model was Slope (Table 5).

The best model for predicting CPUE of all reef species was Area + Trip followed by a model with Period + Trip (Table 2). For assessed reef species CPUE, the same two models were best, though in reverse order. In neither of the top two models did the 95% credible intervals for any predictor variable contain zero in either stage of the model, regardless of whether all reef species CPUE or assessed species CPUE was used (Tables 6, 7). Other models were not competitive.

Qualitative investigations of CPUE revealed interesting trends. When all reef species CPUE was used, a positive temporal trend was evident in the Control area (53.6% increase) but not the SWMPA (9.9% decrease). When only assessed species CPUE was used, the trend for the Control was negative (13.0% decrease) and a positive trend (28.6% increase) for the SWMPA emerged (Figure 3). Overall combined CPUE values (i.e., n fish/n drops across all trips) for assessed species in the SWMPA increased from 0.027 to 0.039 (47.3% increase) from the before to the after period. In the Control area, the same metric decreased from 0.030 to 0.015 (50.4% decrease).

The four Area/Period combinations each had distinct communities. NMDS indicated four groupings with minimal overlap and had an estimated stress value of 0.131, indicating a fair or good fit (Figure 4). NMDS results suggested changes in occurrence from more groupers to more snappers and tilefish in both areas. From the ANOSIM procedure, the R statistic was 0.454 with an associated p value of 0.001, indicating significant separation among the four communities. The results from these multivariate analyses indicate that fish communities have changed without respect to the area (MPA vs. control).

Annual sample sizes for red porgy ages from 2015 to 2020 ranged from 14 to 67 inside the SWMPA and from 8 to 73 inside the control area. For all red porgy from outside the SWMPA, yearly sample sizes ranged from 386 to 546. Mean annual ages ranged from 4.9 to 7.1 inside the SWMPA and from 4.8 to 5.3 in the control area, and were higher in the SWMPA for each of the 6 years examined. Wilcoxon rank sum tests found that ages were significantly higher in the SWMPA than the control in 2016 and 2018 (Figure 5). Age diversity values ranged from 1.7 to 2.0 inside the SWMPA and from 0.9 to 2.0 in the control. Diversity values were higher inside the SWMPA for five out of 6 years examined (Figure 5).

Rudershausen et al. (2010) collected lengths for 26 red porgy in the control area and 41 red porgy in the SWMPA in 2007–2009. Between data from the present study and SERFS data, lengths were available for 89 red porgy in the control area and 181 red porgy in the SWMPA from 2018 to 2020. Size distributions showed an increase in red porgy total lengths in the SWMPA from the before to the after period but no shift was evident in the control area (Figure 6). Our Bayesian linear model revealed that the Area^∗Period interaction was meaningful, as its credible interval did not overlap zero. This result indicates a positive effect of MPA designation on the size structure of red porgy inside the SWMPA.