A marine protected area network does not confer community structure resilience to a marine heatwave across coastal ecosystems

Marine protected areas (MPAs) have gained attention as a conservation tool for enhancing ecosystem resilience to climate change. However, empirical evidence explicitly linking MPAs to enhanced ecological resilience is limited and mixed. To better understand whether MPAs can buffer climate impacts, we tested the resistance and recovery of marine communities to the 2014–2016 Northeast Pacific heatwave in the largest scientifically designed MPA network in the world off the coast of California, United States. The network consists of 124 MPAs (48 no‐take state marine reserves, and 76 partial‐take or special regulation conservation areas) implemented at different times, with full implementation completed in 2012. We compared fish, benthic invertebrate, and macroalgal community structure inside and outside of 13 no‐take MPAs across rocky intertidal, kelp forest, shallow reef, and deep reef nearshore habitats in California's Central Coast region from 2007 to 2020. We also explored whether MPA features, including age, size, depth, proportion rock, historic fishing pressure, habitat diversity and richness, connectivity, and fish biomass response ratios (proxy for ecological performance), conferred climate resilience for kelp forest and rocky intertidal habitats spanning 28 MPAs across the full network. Ecological communities dramatically shifted due to the marine heatwave across all four nearshore habitats, and MPAs did not facilitate habitat‐wide resistance or recovery. Only in protected rocky intertidal habitats did community structure significantly resist marine heatwave impacts. Community shifts were associated with a pronounced decline in the relative proportion of cold water species and an increase in warm water species. MPA features did not explain resistance or recovery to the marine heatwave. Collectively, our findings suggest that MPAs have limited ability to mitigate the impacts of marine heatwaves on community structure. Given that mechanisms of resilience to climate perturbations are complex, there is a clear need to expand assessments of ecosystem‐wide consequences resulting from acute climate‐driven perturbations, and the potential role of regulatory protection in mitigating community structure changes.


| INTRODUC TI ON
Climate change can rapidly reshape the distribution of species and the composition of ecological communities (Pörtner et al., 2023;Smale et al., 2019), imperiling nature's contributions to people. In particular, episodic periods of anomalous ocean warming, hereafter "marine heatwaves," are driving pronounced shifts in species distributions across marine ecosystems (Azzurro & D'Amen, 2022;Olsen et al., 2022), with direct implications for ecological processes and associated human benefits (Cheung et al., 2021;Cinner et al., 2022;Payne et al., 2021;Smale et al., 2019;Smith et al., 2023). While the urgency to plan for adaptation to climate change is clear, as marine heatwaves increase in frequency and severity , pathways to enhance ecosystem resilience are mixed. Therefore, understanding how ecological communities resist, recover from or are transformed by climate perturbations, such as marine heatwaves, represents one of the most pressing challenges for building ecosystem resilience capacity (Mason et al., 2022).
Marine protected areas (MPAs) are an important conservation strategy for preserving biodiversity and nature's contributions to people (Grorud-Colvert et al., 2021;Nowakowski et al., 2023). MPAs may provide network-scale population connectivity that can enhance spillover of individuals and the replenishment of populations in both protected and fished areas (Baetscher et al., 2019;Di Lorenzo et al., 2020;Goetze et al., 2021;Harrison et al., 2012;Williamson et al., 2016). Although most MPAs were initially designed to reduce the effects of overfishing and habitat loss, they are frequently hypothesized to provide long-term protection against climate impacts (Hofmann et al., 2021;Roberts et al., 2017). For example, networks of MPAs may provide refugia as species redistribute in response to climate change, owing to lower anthropogenic stressors and higher population sizes resulting from reduced harvest (Carr et al., 2017;McLeod et al., 2009). MPA features such as habitat diversity, historic fishing pressure, age, and size may also influence the capacity for MPAs to provide ecological resilience (Jacquemont et al., 2022). As such, MPAs and MPA networks are being increasingly highlighted as a key tool for enhancing climate resilience (IUCN-WCPA, 2008;Jacquemont et al., 2022).
Despite the growing number of studies examining MPAs as a tool for mitigating climate impacts, the effectiveness of MPAs (both individually and as networks) for enhancing the resilience (i.e., resistance to and recovery from disturbance, Bates et al., 2019) of marine communities to climate change remains contested (Johnson et al., 2022, but see Jacquemont et al., 2022). Climate change stressors (e.g., ocean acidification, sea level rise, hypoxia, warming) may have effects on populations and communities that occur regardless of regulatory protection Bruno et al., 2019). Marine heatwaves can also reduce connectivity by changing prevailing current patterns-a key design attribute for many MPA networks (Lima et al., 2021). Mixed evidence surrounding the efficacy of MPAs in providing climate resilience may be explained in part by the singlehabitat (e.g., coral reef, seagrass, kelp forest, etc.) focus of many studies, and also in part by the ways in which assemblages of species are partitioned.
Our understanding of whether and how MPAs confer ecological resilience to climate change can be improved by synthesizing the effects of regulatory protection across multiple taxa, habitats, ecosystems, and protection levels. Potential mechanisms producing Arnhold UC Santa Barbara-Conservation International Climate Solutions Collaborative; BiodivERsA; Fondation de France California's Central Coast region from 2007 to 2020. We also explored whether MPA features, including age, size, depth, proportion rock, historic fishing pressure, habitat diversity and richness, connectivity, and fish biomass response ratios (proxy for ecological performance), conferred climate resilience for kelp forest and rocky intertidal habitats spanning 28 MPAs across the full network. Ecological communities dramatically shifted due to the marine heatwave across all four nearshore habitats, and MPAs did not facilitate habitat-wide resistance or recovery. Only in protected rocky intertidal habitats did community structure significantly resist marine heatwave impacts.
Community shifts were associated with a pronounced decline in the relative proportion of cold water species and an increase in warm water species. MPA features did not explain resistance or recovery to the marine heatwave. Collectively, our findings suggest that MPAs have limited ability to mitigate the impacts of marine heatwaves on community structure. Given that mechanisms of resilience to climate perturbations are complex, there is a clear need to expand assessments of ecosystem-wide consequences resulting from acute climate-driven perturbations, and the potential role of regulatory protection in mitigating community structure changes.

K E Y W O R D S
California, climate change, community composition, community structure, marine heatwaves, marine protected area networks, resilience ecological shifts in response to climate change may include altered mortality due to physiological environmental tolerances, changes in species interactions (e.g., competition, predation, disease, facilitation), adult movement across habitats and along the coast, ontogenetic shifts, and changes in recruitment success (Harley et al., 2006). Moreover, the ways in which species responses are evaluated (e.g., trait-based, functional groups, feeding guilds, evolutionary lineages, etc.) can influence detection of climate-driven outcomes. Therefore, long-term monitoring across multiple habitats subjected to similar (or the same) perturbations is needed to thoroughly examine the impacts of climate change on marine ecosystems, to test whether (and which) MPAs confer climate resilience, and to assess the relative resilience among ecosystems and taxa. Such cross-ecosystem studies are rare, especially those that include monitoring across networks of MPAs before, during, and after extreme climate change-driven perturbations.
In 1999, California passed the Marine Life Protection Act (MLPA), which expanded its system of MPAs to function as a coherent ecological network and to address six goals aimed at conservation, fisheries, and other human benefits (Gleason et al., 2013;Marine Life Protection Act, 1999). Guided by these goals, California established a network of 124 MPAs (48 no-take state marine reserves, 76 partial-take or special regulation conservation areas) distributed along the state's entire 1300-km coastline that protects 16% of state waters. The network protects hard-and soft-bottom habitats ranging in depth from the intertidal to depths of 1000 m. However, with few exceptions, ecological monitoring studies focused on hard-bottom habitats. Leading up to and following MPA implementation, an extensive ecological monitoring effort of these habitats began to support adaptive management of the network (Botsford et al., 2014). During the course of California's MPA monitoring, a major marine heatwave occurred that was the consequence of two environmental anomalies: a 2014-2015 warming event known as "the Blob," and a major El Niño event in 2015-2016 (Bond et al., 2015;Di Lorenzo & Mantua, 2016;Gentemann et al., 2017). This pronounced climate perturbation was the largest marine heatwave on record in California (Laufkötter et al., 2020).
This study leveraged over a decade of monitoring across the rocky intertidal, kelp forest, shallow reef, and deep reef habitats to evaluate whether a network of MPAs confers community structure resilience to climate disturbances. Specifically, we tested the following hypotheses: (1) the 2014-2016 marine heatwave event resulted in the reorganization of community structure across multiple nearshore habitats and taxa, (2) MPAs enhanced the ability for ecological communities to resist and recover from the impacts of climate perturbations, and (3) species traits (thermal affinities) explained differential responses to the marine heatwave in relation to regulatory protection (MPAs). After testing these hypotheses, we explored whether MPA features (size, age, historic fishing pressure, habitat diversity, and connectivity) enhanced ecological community resilience across habitats to inform the design of MPA networks that are resilient to climate change.

| MPA sampling and data collection
As part of the State of California's long-term MPA evaluation and monitoring program, several habitat-specific research groups conduct annual surveys designed to monitor ecological changes over time inside MPAs and at areas of comparable habitat outside MPAs ( Figure 1a; Figure S1). Conceptually, while MPAs could be considered the "treatment" or "control" in our analyses (i.e., fishing is considered either the experiment or the control), we take a social-ecological system perspective and consider human activities, including fishing, to be integral system components (Kandel et al., 2022). Thus, in our analytical framework, MPAs are considered the experimental treatment where fishing is removed (e.g., Pacoureau et al., 2023).
Our analyses focused on four habitats that have extensive longterm monitoring data: rocky intertidal, kelp forest, shallow reef, and deep reef (Figure 1a; Figure S1). Across the four habitats included in our analyses, three organismal groups were sampled: macroalgae (rocky intertidal and kelp forest), conspicuous mobile and sessile invertebrates (rocky intertidal and kelp forest), and fishes (kelp forest, shallow reef, and deep reef). We focused our community structure analyses on 13 MPAs located along the temperate Central Coast of California ( Figure 1a) because this area was the most comprehensively sampled region of the MPA network and had sufficient premarine heatwave data to evaluate baseline community structure (Table S1; Figure S1). To evaluate MPA features as potential drivers of ecological resilience, we included 28 MPAs across all of California for two of the habitats (rocky intertidal and kelp forest), where more extensive spatial and temporal coverage exists (see Section 2.6).
Two general types of data were used in our analyses: species counts (kelp forest, shallow reef, deep reef) and the proportional cover of invertebrates and macroalgae (rocky intertidal and kelp F I G U R E 1 The (a) coverage of ecological monitoring in marine protected areas (MPAs) in California's Central Coast region and (b-f) exposure of these MPAs to five indicators of environmental conditions. In (a), points indicate which habitats were monitored both inside and outside of MPAs (i.e., points indicate data availability, not the location of sampling sites). Red polygons indicate no-take MPAs (n = 12, called State Marine Reserves), and blue polygons indicate partial-take MPAs (n = 10, called State Marine Conservation Areas; see Table S1 for list of partial-take MPAs that are de facto no-take MPAs by habitat). The dark horizontal lines delineate the Central Coast region, and the thin gray line indicates state waters (three nautical miles offshore and all of Monterey Bay). Environmental indicators include (b) air temperature (AT); (c) sea surface temperature (SST); (d) sea bottom temperature (SBT); (e) the multivariate oceanographic climate index (MOCI); and (f) the biologically effective upwelling index (BEUTI). Lines indicate the median and shading indicates the 95% confidence interval. The 2014-2016 marine heatwave (MHW) is indicated by the vertical red rectangle.
forest invertebrates and algae). While biomass is often used as a common measurable response across taxa (Duffy et al., 2017), it is not suitable for our analyses because of the large number of taxa for which it is difficult to accurately measure biomass (e.g., macroalgae and invertebrates). Therefore, we elected to use a taxonomic abundance-based approach to compare changes in community structure. All analyses were conducted using published data for each habitat Cieri et al., 2022;Malone et al., 2022;MARINe et al., 2022).

| Changes in community structure associated with the marine heatwave
We used non-metric multidimensional scaling (nMDS) to visualize community structure before (2007-2013), during (2014-2016), and after (2017-2020) the 2014-2016 marine heatwave. Prior to ordination, a base similarity matrix was constructed for each habitat type using Bray-Curtis dissimilarity on Hellinger-transformed counts or on percent cover of species. The Hellinger transformation converts absolute counts to the square root of proportional counts, which reduces the disproportionate contribution of highly abundant species (Legendre & Gallagher, 2001). We elected to use the transformation to relative species abundance to improve the robustness, comparability, and ecological interpretation of our community structure analyses. Each site sampled inside or outside of a MPA surveyed in a single year represents a single unit of replication. To visualize the state of each ecological community, we plotted centroids and 95% confidence ellipses that represent the generalized position of the community in ordinated two-dimensional nMDS space before, during, and after the marine heatwave event.
We used a pairwise permutational analysis of variance (PERMANOVA) to test three hypotheses surrounding community structure changes ( Figure 2a). First, if MPAs do not mitigate marine heatwave impacts on community structure, then communities inside MPAs should respond similarly (i.e., comparable change in multivariate distance) to those outside. Second, if MPAs confer resistance, then communities inside MPAs should remain unchanged while those outside shift. Finally, if MPAs enhance recovery, then communities inside MPAs should shift during the marine heatwave, but return to or move back in the direction of their previous state following the disturbance. To test these hypotheses, a single PERMANOVA was performed for each habitat and site type (inside or outside a MPA) using pairwise comparisons between marine heatwave periods (before vs. during, before vs. after) using the pairwiseAdonis wrapper function (Martinez Arbizu, 2020) in the vegan package (Oksanen et al., 2022) in R (R Core Team, 2021).

| Effect of MPAs on community structure resistance and recovery
We used a multivariate distance-based approach to test whether MPAs conferred ecological community resistance to or recovery from the marine heatwave. We refer to this measurement as "multivariate distance" to distinguish it from "geographic distance," which is often a measure of climate change impacts on marine species distributions. We define resistance as community structure (i.e., relative abundance of species) that remained unchanged during the marine heatwave (low change in multivariate distance-based centroids), and recovery as a community that returned to a similar structure postheatwave. Resistance was evaluated by calculating the multivariate vector distance (in high-dimensional space using the Bray-Curtis dissimilarity matrix) of the centroid of the ecological community between the periods before (2007)(2008)(2009)(2010)(2011)(2012)(2013) and during (2014)(2015)(2016) the marine heatwave (e.g., high resistance is indicated by a smaller change in multivariate distance-based centroids).
Recovery was evaluated by calculating the multivariate distance between the before (2007-2013) and after period (2017)(2018)(2019)(2020) centroids. We calculated the change in multivariate distance between centroids using the betadisper function in the vegan package (Oksanen et al., 2022). The betadisper function returns the principal coordinates of centroids, which we used to calculate multivariate distance after ensuring positive-definite eigenvalues. Finally, statistical significance was evaluated using a PERMANOVA on community structure resistance (before vs. during) and recovery (before vs. after).
To explore whether the timing of community shifts coincided with temporal changes in oceanographic variables, we used a Granger causality test on Bray-Curtis dissimilarity. For this analysis, annual dissimilarity was calculated for each site (inside or outside an MPA) and habitat relative to 2007. Dissimilarity was then offset (lagged) against oceanographic variables (see Section 2.5) in 1-year increments for a maximum of 3 years (maximum lag based on length of the time series). We used 2007 as the baseline year because it preceded the marine heatwave and because we were interested in lag effects specifically related to the marine heatwave, rather than gradual environmental changes over time. However, 2008 was used as the baseline year for the deep reef habitat because surveys were not extensively conducted in 2007. It is important to note that the Granger test only examines lagged effects between the time series of community change and the environmental variables, it does not examine correlations, which are explored using other diagnostics described in Section 2.5.

| Species associated with community changes
To evaluate which species were associated with community change, we used two approaches. First, we used a similarity percentage analysis (SIMPER) to decompose community structure using Bray-Curtis dissimilarity and to examine the percent contribution of each species to the before versus after marine heatwave communities.
However, because SIMPER is known to confound the mean between group variation and dispersion (Warton et al., 2012), we crossvalidated the SIMPER output using the mvabund package (Wang et al., 2022) in R. The mvabund package uses fitted generalized linear models to account for nonlinear mean-to-variance relationships of each species. Both of these analyses were conducted on the absolute abundance of species (rather than relative abundance) because they are not sensitive to common or rare taxa. This combination of approaches allowed us to determine the contribution of individual F I G U R E 2 (a) Potential and (b) observed shifts in community structure from 2007 to 2020. (a) Illustrates a typology of potential community structure shifts: (i) "No impact," where community structure does not shift, and "impact", where community structure significantly shifts. Non-overlapping ellipses indicate significant differences in community structure. In this example, pre-heatwave community structure was significantly different inside (solid lines ellipses around points) versus outside (dashed lines around triangles) MPAs; (ii) "No MPA benefit," where communities inside MPAs respond similar to those outside; (iii) "MPA Resistance," where communities inside MPAs remain stable while those outside shift; and (iv) "Recovery," where communities inside MPAs shift during the heatwave, but return to their previous state following the disturbance while communities outside of MPAs do not return. (b) Shows observed shifts in community structure before (2007-2013), during (2014-2016), and after (2017-2020) the marine heatwave inside (circles) and outside (triangles) MPAs in the Central Coast region of California. Each subpanel represents a habitat and each point depicts the centroid position with 95% confidence ellipses. Also included are vectors for each environmental anomaly: air temperature (AT, rocky intertidal only), sea surface temperature (SST), sea bottom temperature (SBT), the multivariate ocean climate index (MOCI), and the biologically effective upwelling transport index (BEUTI). The trajectory of each vector reflects its correlation with community structure (Table 2). Therefore, indicators that are highly correlated with changes in community structure are aligned with the centroids (points). MPA, marine protected area; nMDS, non-metric multidimensional scaling. taxa to observed community structure differences. Additionally, because the PERMANOVA analyses revealed that community change was similar inside and outside of MPAs, we analyzed compositional changes overall rather than within each site by protection status.
Finally, we constructed an analysis of deviance table using the multivariate generalized linear model fits to test whether species structure changed as a result of the marine heatwave (before vs. after).
Test statistics and p-values for each species were generated using the PIT-trap (probability integral transform) resampling method.

| Environmental correlates and species traits
We explored temperature and oceanographic conditions before, during, and after the marine heatwave to evaluate whether community structure shifts were explained by environmental changes indicator of several oceanographic and atmospheric conditions (García-Reyes & Sydeman, 2017) calculated at the regional level (Central California). We selected these environmental indicators and associated products because they provide biologically meaningful climatology at the best available spatial resolution for our study area.
To process the environmental data, we first calculated the monthly mean AT, SST, SBT, BEUTI, and quarterly MOCI values at each site. We then calculated monthly anomalies for AT, SST, SBT, and BEUTI as the difference between the observed monthly mean

| Drivers of resistance and recovery
The ecological and design context of California's MPA network is well documented, allowing us to explore whether certain MPA features conferred community structure resilience to the marine heatwave. These analyses were spatially expanded to include a subset of MPAs covering the entire network from northern to southern California, thereby providing a greater scope of MPA feature diversity ( Figure S3). However, only the kelp forest and rocky intertidal habitats contained sufficient pre-heatwave monitoring data with consistent annual surveys in all bioregions of the state to appropriately examine drivers of ecological resilience among MPAs at the network level (Figure 1a; Figure S1). MPA features included habitat richness and habitat diversity (calculated using Shannon-Wiener on depth-stratified hard and soft bottom, proportion of rock, and the extent of kelp forest canopy, rocky intertidal, sandy beach, coastal marsh, tidal flats, and armored shore), historic fishing pressure, MPA age and size, connectivity, and kelp forest fish biomass response ratios (as a proxy for MPA ecological performance, see Supplementary Methods, Appendix S1).
We developed a two-stage multivariate model to examine the effect of MPA features on community resistance and recovery. First, we used a PERMANOVA to identify MPAs that significantly resisted or recovered from the marine heatwave (building on the methods outlined in Section 2.3). We then used logistic regression to evaluate the probability that specific MPA features enhanced resistance or resilience. For the logistic regression, any MPA community that significantly resisted or recovered from the marine heatwave was assigned a target level of "1," and any MPA community that did not resist or recover from the marine heatwave a "0."

| Community structure resistance and recovery
Ecological community change was widespread across all habitats, and three habitats (kelp forest, shallow reef, deep reef) showed no clear differences inside and outside of MPAs in the magnitude of community change (Figure 2b). Prior to the marine heatwave, community structure was similar inside and outside of MPAs in all habitats except the rocky intertidal (indicated by the non-overlapping ellipses in Figure 2b). However, ecological community structure dramatically shifted as a result of the marine heatwave event across all measured nearshore habitats (n = 4), regardless of whether the communities (n = 5) were inside or outside of MPAs (Figure 2b; Figure S8; Table 1).
Community structure changes coincided with oceanographic conditions associated with the marine heatwave (Figure 1; Tables 1 and 2). At the onset of the marine heatwave in 2014, SST was anomalously warm by as much as 2°C, and sea bottom temperature was also 1°C above the baseline average. MOCI and BEUTI experienced precipitous changes during the marine heatwave, reflecting reduced upwelling and productivity, and these anomalies persisted until at least mid-2016 (Figure 1e,f). Although each oceanographic anomaly aligned well with community shifts in nMDS space (Figure 2b), significant correlations were habitat-specific (Table 2). In the rocky intertidal, warmer AT was the strongest environmental correlate with community structure changes, although warmer SST and declines in upwelling (BEUTI) were also correlated with community change.
For shallow reef, kelp forest fishes, and kelp forest invertebrates and algae, SST, SBT, and MOCI were significantly associated with community change, such that communities shifted with increases in surface and bottom water temperature and positive values of the MOCI index. In the deep reef habitat (fishes only), the only significant oceanographic correlate was BEUTI (Table 2), where the community shifted in response to reduced upwelling. Finally, there were no significant temporal lags determined except for one interaction between the shallow reef and BEUTI (lag = 3 years for sites inside and outside of MPAs, Table S4).
For all habitats except the rocky intertidal, MPAs did not impart increased resistance to or recovery from marine heatwave-driven community changes compared to sites outside of MPAs (Figure 3).
The rocky intertidal and deep reef were the only two habitats that resisted the marine heatwave (based on nonsignificant PERMANOVA result, Figure 3a). The rocky intertidal was the only habitat where community structure inside MPAs was not significantly different either during or post-marine heatwave (Figure 3b). Kelp forest invertebrates, algae, and fishes all substantially shifted during the marine heatwave, and these changes persisted for several years later with no apparent recovery. The shallow reef habitat responded similarly with pronounced shifts in community structure during the marine heatwave, although the trajectory of the shallow reef habitat started to move toward the pre-heatwave state beginning in the year 2018, with the biggest shifts toward recovery occurring inside MPAs ( Figures S5 and S6). In the deep reef habitat, year-to-year changes TA B L E 1 Results from a series of pairwise permutational analysis of variance (PERMANOVA) tests on community resistance (beforeduring) and recovery (before-after). PERMANOVAs were performed separately for each habitat and calculated on Bray-Curtis dissimilarity matrices. were more variable, resulting in a larger multivariate distance relative to the other habitats, but this community change was not significantly different from the baseline year (Figure 3a; Figure S9).
However, recovery was more variable than resistance. The shallow reef habitat exhibited relatively greater recovery (i.e., less shift in multivariate distance compared to the pre-heatwave community state) than the kelp forest and deep reef habitats.

| Species responses
The multivariate analyses revealed several species that explained differences between the pre-and post-heatwave periods. Among the three habitat types with monitoring of fish species (kelp forest, shallow reef, deep reef), the blue and deacon rockfish complex (Sebastes mystinus and S. diaconus) was positively correlated with the post-heatwave period (

| Community structure and thermal traits
The relative proportional representation of thermal affinities significantly changed during the marine heatwave for kelp forest invertebrates and algae, kelp forest fishes, and shallow reef fishes (p < 0.001, F 6,708 = 62.24; p < .05, F 6,700 = 2.18; p < .001, F 4,315 = 37.41, respectively). Cold-temperate species significantly declined during the marine heatwave for these habitats (Figure 4a). During the marine heatwave, there was a slight increase in cosmopolitan species for the rocky intertidal and kelp forest (invertebrates and algae), and an even more pronounced increase in warm-temperate and subtropical fish species in the kelp forest, shallow reef, and deep reef habitats. Importantly, these changes in community composition that occurred during the marine heatwave (2014-2016) persisted into the following years (2017 and beyond), which partially explains the lack of observed recovery to the pre-heatwave community structure.
Variation in species abundance was explained by significant interactions between thermal affinities and oceanographic variables for all habitats from the fourth-corner model (Figure 4b; Figure S7).

| MPA features and ecological stability
Community structure responses were highly variable across the state-level network for MPAs. In general, rocky intertidal communities showed higher resistance than kelp forest communities, as indicated by smaller shifts in multivariate distance (Figure 5c; smaller points). While shifts in community structure were generally larger in the kelp forest (Figure 5c; larger points), in some cases, the magnitude of shifts inside MPAs was less relative to shifts outside of MPAs (Figure 5a). In general, MPA protection fostered more resistance than recovery (Figure 5c; i.e., more blue points in the resistance column and more red points in the recovery column), although these results were not statistically significant. For example, kelp forest communities in the Point Buchon, South Point, and Harris Point MPAs resisted and moved toward recovery from the heatwave more than communities outside these MPAs ( Figure 5; blue points).
However, kelp forest communities in the Scorpion, Anacapa Island, and Campus Point MPAs were less successful at resisting shifts than communities outside these MPAs (Figure 5c, left; red points).
Similarly, kelp forest communities in the Point Lobos, Santa Barbara, Farnsworth, and Blue Cavern MPAs were less successful at recovering from shifts than communities outside these MPAs (Figure 5c, right; red points). None of the evaluated MPA features (age, size, historic fishing pressure, habitat diversity, connectivity, biomass response ratios) were statistically significant drivers of this variability ( Figure S11).

| DISCUSS ION
Our findings suggest that the 2014-2016 Northeast Pacific marine heatwave impacted ecological community structure across four nearshore habitats, and that MPAs did not confer widespread resistance or recovery. MPAs have gained increased attention as a conservation strategy for mitigating the effects of climate change Wilson et al., 2018), but evidence of their efficacy in providing ecosystem-wide resilience to climate disturbances remains mixed Jacquemont et al., 2022;Roberts et al., 2017). Critically, international efforts to conserve 30% of marine habitats via MPA implementation by the year 2030 (Day et al., 2012) highlight the need for planning that considers the effect of large climate-driven perturbations on local ecological processes.
Therefore, while MPAs may be implemented with specific conservation targets, our study suggests that extreme climate perturbations, such as marine heatwaves, can overwhelm intended climate benefits of MPAs over the short-term (<5 years following a marine heatwave).
These results highlight the ecosystem-wide consequences of acute climate-driven perturbations despite regulatory protection.

F I G U R E 3
Community resistance [(a) before to during marine heatwave comparison] and recovery [(b) before to after marine heatwave comparison] as measured by multivariate distance between centroids. Each point depicts the multivariate distance between the preheatwave centroid (before) and the during-heatwave or the post-heatwave centroid inside (red) and outside (blue) marine protected areas. Therefore, higher values indicate less resistance or recovery. Error bars depict the pooled standard error between centroids and the asterisks denote significant differences in community structure between marine heatwave periods, as derived from a pairwise permutational analysis of variance (PERMANOVA). Points without asterisks indicate resistance or recovery (as measured by nonsignificant PERMANOVA, Table 1). Note that some communities have inherently greater interannual variability, potentially resulting from variation in sampling methodology, which can lead to larger mean multivariate distances that are not significantly different ( Figure S9).
Several ecological processes may explain the pronounced shifts

F I G U R E 5
Conceptual schematics illustrating scenarios where (a) community structure shifts less inside an MPA than outside, thereby reducing community structure shifts, and (b) community structure shifts more inside an MPA than outside, thereby exacerbating community structure shifts relative to the outside sites. (c) Shows resistance (before-to-during marine heatwave comparison) and recovery (before-toafter marine heatwave comparison) of community structure to the 2014-2016 marine heatwave by MPA and habitat type along a latitudinal gradient. Point size indicates the Bray-Curtis multivariate distance for resistance or recovery. Therefore, smaller points indicate MPAs that exhibited greater resistance or recovery (i.e., communities inside the MPA shifted smaller multivariate distances). Point color indicates the magnitude of the shift of community structure inside a given MPA relative to its paired outside site. Red shades indicate MPAs where the change in multivariate distance was greater than the paired outside site (i.e., shift exacerbated inside MPAs relative to outside). Blue shades indicate MPAs where the change in multivariate distance was greater in the paired outside site than inside the MPA (i.e., community structure changes were less in the MPA). Note that the rocky intertidal does not use a paired inside-outside sampling design (outside site selection is described in Supplementary Methods, Appendix S1), and therefore, the grey points indicate MPAs without the necessary data from adjacent outside sites. Missing points indicate insufficient data to evaluate changes. MPA, marine protected area; SMCA, State Marine Conservation; SMR, State Marine Reserve (no-take MPA). study found that the relative proportion of warm-temperate and subtropical species increased during and after the marine heatwave.
The increased proportional representation of species with warm water thermal affinities may be explained by adult movement into the study area, changes in recruitment patterns, or by a decline in the relative abundance of cold-temperate species (Fredston et al., 2021;Sanford et al., 2019;Walker et al., 2020). For example, the large increase in blue rockfish in kelp forest, shallow reef, and deep reef habitats was associated with strong recruitment and pelagic youngof-the-year abundance of rockfishes in midwater trawl surveys at the start of the marine heatwave period (Field et al., 2021). Large increases in the subtropical wrasse, senorita (Oxyjulis californica), were also observed in multiple habitats, contributing to community change. Second, during the marine heatwave, a large-scale outbreak of herbivorous sea urchins occurred throughout the study region that coincided with precipitous declines in kelp (and other macroalgae) and the loss of an important benthic mesopredator, the sunflower sea star (Pycnopodia helianthoides), due to a large-scale marine disease outbreak (Harvell et al., 2019;Smith et al., 2021). The primary regulation associated with MPAs in California involves a restriction of fishing activities; thus, we were interested in understanding whether ecological communities inside no-take MPAs are more resilient to marine heatwaves than our control social-ecological system. We consider human activities, such as fishing, as integral parts of the social-ecological system, since fishing occurred decades before MPA establishment and the marine heatwave, with the restriction of harvest through regulatory protection as the primary conservation mechanism. However, there are multiple viewpoints pertaining to interactions between humans and the oceans, and similarly, there are many approaches to evaluating the effects of fishing and regulatory protection. Alternative frameworks might consider areas inside MPAs as "reference" locations for places where harvest occurs (Costello, 2014). While we believe human and natural systems are now so integrated in coastal areas that there are virtually no pristine ecosystems, the framing of considering MPAs as controls and fished systems as experimental treatments would not impact our results.
It is not surprising that community shifts occurred inside and outside MPAs simultaneously for most habitats, since many of the species exhibiting the biggest changes are not directly targeted by fishing activities (e.g., invertebrates and algae). In cases where MPAs do not confer ecological resilience to marine heatwaves, it is important to note that the capacity of MPAs to harbor higher biomass than unprotected sites can still be maintained (Frid et al., 2023). Other studies have documented increases in targeted species biomass inside of MPAs within this system, which is a focal MPA conservation objective (Caselle et al., 2015;Hamilton et al., 2010;Ziegler et al., 2022). Similar responses to the marine heatwave in southern California fish communities were reported by Freedman et al.
(2020), who found that non-targeted subtropical species were most responsible for community reorganization. Ziegler et al. (2023) also reported that fish species diversity recovered more quickly after the heatwave inside of MPAs.
Although the anomalously warm environmental conditions subsided after 2016, changes in community structure persisted, especially in the kelp forest, shallow reef, and deep reef habitats. These results are particularly interesting in the context of ecosystem stability and transition dynamics. Given the highly variable recruitment of the dominant fishes in the study region (Field et al., 2021;Schroeder et al., 2019), it is possible that sufficient time has not elapsed for the communities to return to a pre-perturbed state. Lagged effects due to ontogenetic shifts and methodological sampling that disproportionately select for adults may also contribute to the observed lack of recovery. For example, declines in vermilion rockfish (S. miniatus) in the shallow reef habitat coincided with increases in that species in the deep reef habitat, likely reflecting ontogenetic movements from shallow to deep habitats (Love, 2011;Love et al., 2002). In addition, the shallow reef habitat monitoring program used hook and line sampling, which disproportionately selects individuals that are older than 2-3 years, which is why lagged effects of oceanography were more likely to be detected in this habitat (Ziegler et al., 2023).
Moreover, many species in the California Current are slow-growing, late to mature, long-lived, and relatively sedentary with small home ranges (Love, 2011;Love et al., 2002). These life-history traits may create inertia to resistance or recovery from change. However, the persistence of novel community structure configurations observed in this study, despite the return of pre-heatwave environmental conditions, highlights the need for further assessments of community transition dynamics over ecologically meaningful timescales.
Our finding of no effect of MPA features on dampening marine heatwave-driven responses in the rocky intertidal and kelp forest habitats is particularly interesting given that habitat mosaics, historic fishing effort, and connectivity have been linked to MPA conservation capacity (Bastari et al., 2016). The widespread effects of marine heatwaves and other climate perturbations may overwhelm local stressors (e.g., fishing, marine pollution) that MPAs were initially successful at mitigating . Indeed, even biodiversity-rich spots or deep water communities are susceptible to the impacts of climate change (Emblemsvåg et al., 2022;Kocsis et al., 2021). Our results highlight the need for management frameworks that include climate adaptation strategies, the protection of range-shift corridors, and consideration of the compounding effects of both fishing and marine heatwaves on resilience (Burrows et al., 2014). More investigation is needed to better understand the drivers that can help MPAs to maintain higher fish biomass than unprotected sites under climate change (Frid et al., 2023).
Increasing our understanding of the pathways through which marine heatwaves restructure ecological communities is central to developing adaptive management solutions. In our study, we used a taxonomic (i.e., species) abundance-based approach to evaluate changes in community structure, but functional diversity may also be impacted by warming events and should be explored (McLean et al., 2019;Murgier et al., 2021). Prolonged warming events may differentially impact groups and guilds of species with similar traits, such as dispersal ability, thermal tolerance, metabolic rate, and mobility Harvey et al., 2021). Moreover, to understand the ability of MPA networks to resist and recover from future marine heatwaves, it is critical to have sufficient monitoring across multiple habitats, taxa, and regions. Because pre-perturbation data are essential for establishing baselines and comparatively evaluating ecological responses, consistent monitoring is fundamental to capture the effect of marine heatwave events. Although the timing of marine heatwaves is unpredictable, evidence suggests that the frequency and magnitude of abrupt warming events are expected to increase (Frölicher et al., 2018;Holbrook et al., 2020).

ACK N OWLED G M ENTS
We thank the hundreds of individuals who contributed to data collection for each of the long-term monitoring programs included in this study. We also thank R. Brooks

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors have no conflicts of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly