Cataclysmic Disturbances to an Intertidal Ecosystem: Loss of Ecological Infrastructure Slows Recovery of Biogenic Habitats and Diversity

Understanding the resilience and recovery processes of coastal marine ecosystems is of increasing importance in the face of increasing disturbances and stressors. Large-scale, catastrophic events can re-set the structure and functioning of ecosystems, and potentially lead to different stable states. Such an event occurred in south-eastern New Zealand when a Mw 7.8 earthquake lifted the coastline by up to 6 m. This caused widespread mortality of intertidal algal and invertebrate communities over 130 km of coast. This study involved structured and detailed sampling of three intertidal zones at 16 sites nested into four degree of uplift (none, 0.4–1, 1.5–2.5, and 4.5–6 m). Recovery of large brown algal assemblages, the canopy species of which were almost entirely fucoids, were devastated by the uplift, and recovery after 4 years was generally poor except at sites with < 1 m of uplift. The physical infrastructural changes to reefs were severe, with intertidal emersion temperatures frequently above 35°C and up to 50°C, which was lethal to remnant populations and recruiting algae. Erosion of the reefs composed of soft sedimentary rocks was severe. Shifting sand and gravel covered some lower reef areas during storms, and the nearshore light environment was frequently below compensation points for algal production, especially for the largest fucoid Durvillaea antarctica/poha. Low uplift sites recovered much of their pre-earthquake assemblages, but only in the low tidal zone. The mid and high tidal zones of all uplifted sites remained depauperate. Fucoids recruited well in the low zone of low uplift sites but then were affected by a severe heat wave a year after the earthquake that reduced their cover. This was followed by a great increase in fleshy red algae, which then precluded recruitment of large brown algae. The interactions of species’ life histories and the altered physical and ecological infrastructure on which they rely are instructive for attempts to lessen manageable stressors in coastal environments and help future-proof against the effects of compounded impacts.


INTRODUCTION
There has been renewed interest in large and usually infrequent disturbances. As the pace of change in climatic extremes has increased, there is an increasingly urgent need to better understand processes underlying the persistence, resistance, resilience, adaptive capacity and recovery from events such as catastrophic fires, heat waves, hurricanes, and floods across all of earth's ecosystems (Lubchenco and Karl, 2012). Such understanding is a fundamental requirement for the design of effective management (Dale et al., 1998;Chambers et al., 2019). The basic concept involves attention to the vulnerabilities and resilience of ecosystems and communities. Holling's (1973) initial definition of resilience was a measure of the ability of systems to absorb changes of state and driving variables, and yet still persist, and stability as the ability of a system to return to an equilibrium state after a temporary disturbance. Many definitions and the implications of component processes and thresholds have been a focus of subsequent discussions (e.g., Walker et al., 2004;Brand and Jax, 2007). The severity or intensity, frequency or return time, and scale of disturbances have been repeatedly highlighted as features that underpin the robustness of ecosystems and their ability to resist and recover from disturbance (e.g., Levin and Lubchenco, 2008). Large disturbances can re-set the "successional clock, " leaving a residual assemblage that provides a legacy on which subsequent patterns build (Paine et al., 1998). Falk et al. (2019) state that resilience responses are emergent properties from the component processes of persistence, recovery and re-organization. Recovery is influenced by the spatial and temporal scale of initial impacts, as well as the interactions and requirements of key species in the post-disturbance environment. Depending on all of these factors, an ecosystem may or may not return to its former state.
Many recovery outcomes depend on the extent to which the fundamental ecological "infrastructure" has been altered by a disturbance. The concept of "ecosystems as infrastructure" has been discussed since the 1980s and has many definitions and uses involving the elements of inter-related systems providing goods and services to humans (Fulmer, 2009;da Silva and Wheeler, 2017). The concept can be extended, however, to include inter-relationships between species, and the provision of conditions conducive to the survival and persistence of the ecosystems themselves. These "goods and services" of their own making underpin the health and resilience of many characteristic ecosystems, particularly those that have become adapted to survive and thrive in harsh environments.
In the case of intertidal marine systems, much of the biophysical infrastructure on which species and communities build is affected by a wide range of stressors operating over many spatial and temporal scales. The litany of stressors is long and includes eutrophication, marine heat waves, overfishing, invasive species, storms and wave events, coastal development, and sediment from intensive land use, among others (e.g., Schiel, 2009), and additional step-change events (e.g., Orchard et al., 2021). Each of these is known to affect coastal ecosystems through impacts on vulnerable species, ecological structure, and diversity, thereby altering ecological functions and services. Although highly context-dependent because of varying species' life history traits, life spans and ecological niches, natural recovery is compromised by the degree of change in the physical conditions on which they rely and to which they have adapted (e.g., Bekkby et al., 2020). On nearshore rocky reefs, such changes may involve increased air or sea temperatures (Schiel et al., 2004;Cavanaugh et al., 2019;Smale et al., 2019;Thomsen et al., 2019), smothering by sediments (Airoldi, 2003;Schiel, 2006), a compromised light environment in the water column (Tait et al., 2021), human access and direct impacts on reef communities (Povey and Keough, 1991;Schiel and Taylor, 1999;Van De Werfhorst and Pearse, 2007), and altered wave forces (Denny, 1985;Gaylord and Denny, 1997;Schiel et al., 2016), all of which can affect the fundamental requirements of benthic species to attach, survive and thrive.
These issues came to the fore when a large earthquake struck the northeast coast of the South Island of New Zealand (NZ) in November 2016. The Mw 7.8 Kaikôura earthquake event was one of the most complex ever recorded Holden et al., 2017;Shi et al., 2017;Xu et al., 2018). Although originating inland, the stress release activated a large number of faults in a northeasterly direction affecting extensive areas of land and sea Hamling et al., 2017;Gusman et al., 2018). Vertical displacement affected over 130 km of coastline in a highly variable manner but predominantly in the direction of uplift Orchard et al., 2021). Impacts on anthropogenic infrastructure included severe damage to road and rail networks that run close to the sea in this area (Kaiser et al., 2017). Major changes to relative sea-levels, associated with displacement, caused long-lasting effects in comparison to the tidal range of c.2 m. Broad-scale land and seascape changes included the generation of over 170 ha of new terrestrial land and extensive remodeling of intertidal ecosystems (Orchard et al., 2021). Environmental impacts included widespread mortality of algae, invertebrates and fish along virtually the entire coastline, which effectively reset the nearshore ecosystem . Some remnant populations remained in areas of lesser uplift but the natural stage was mostly a blank slate for the ecological communities to re-assemble (Orchard et al., 2021). Ecological resilience was greatly tested, with "recovery" to the former state being far from certain because of the loss of connectivity between remaining patches of key habitat-formers such as large algal species. Previous small-scale experiments in which dominant algal canopies had been removed showed that it took up to 8 years for canopies to re-develop and associated communities to re-assemble, even though clearances of < 10 m 2 were surrounded by reproductively active adults (Lilley and Schiel, 2006;Schiel and Lilley, 2011).
Initial post-earthquake surveys showed that some species of large brown algae, almost all of them fucoids, were functionally extinct along many areas of the uplifted coast, so there were few local sources of propagules to re-colonize affected areas Thomsen et al., 2019). It was anticipated, therefore, that re-establishment of algal assemblages would take many years, and that there would likely be a successional sequence as canopies re-formed and understory species that relied on facilitative effects of canopy cover eventually became established (Schiel and Lilley, 2007). There was also ongoing disturbance from post-earthquake shifting of gravels and sand that buried some reefs, while erosion effects were prominent at others (Orchard et al., 2021). Invertebrate populations suffered high mortality with most limpets, and turbinid and trochid gastropods dying relatively soon after the earthquake. Of particular importance economically was the great mortality of NZ black-foot abalone, Haliotis iris (Gerrity et al., 2020). Tens of thousands of individuals were propelled upwards beyond the tidal influence where they suffered severe heat stress and died (picture in Schiel et al., 2019). Their recruitment habitats among small rocks in the lowest intertidal and upper subtidal zone were also uplifted along much of the coast, causing concerns about future recruitment rates.
The earthquake, therefore, provided a novel opportunity to test recovery dynamics of a severe and rare large-scale disturbance on a coastal marine system. Other similar studies, primarily from earthquake-prone Chile (e.g., Castilla, 1988;Castilla et al., 2010;Jaramillo et al., 2012;Ortega et al., 2014) and Japan (Sato and Chiba, 2016;Muraoka et al., 2017) have shown multi-year effects on the recovery of both soft shore and rocky reef communities.
The initial hypothesis of our study was that the intertidal zone would eventually shift downward to re-establish on newly available boulders and reefs pushed up from the subtidal zone. We set out to test this with intensive, field-based quantitative surveys done annually for 4 years after the earthquake. To relate recovery dynamics and trajectories to physical changes, we monitored temperatures at many sites, gauged erosion and breakup of the soft sedimentary rocks and substrate accretion effects, and monitored the nearshore light environment. We evaluated changes across tidal zones relative to the initial community structure, which was known from surveys and prior data. We also tested to what extent recovery was related to the degree of uplift, as a gauge of the degree of initial disturbance and relative sealevel change. The results facilitate an evaluation of the recovery dynamics and resilience of this complex reef system.

Sites and Design
The Mw 7.8 earthquake struck just after midnight, at low tide on 14 November 2016, near the start of austral summer. The coastal zone was lifted by up to c. 6 m around Waipapa Bay, over 2 m near Cape Campbell in the north, < 1 m on Kaikôura Peninsula, and around 1.7 m at Omihi further south with a considerable degree of variation between sites in localized areas (Orchard et al., 2021; Figure 1). Coastal cliffs fell over the main coastal highway in many places, blocking access from the north and south, and roads remained closed for over a year. The coast is sparsely populated except around the town of Kaikôura, which is a major tourism destination, and access to coastal sites was difficult. There were few detailed data for intertidal and subtidal communities along this coast except for sites to the north and south, which have been surveyed annually for > 20 y (Schiel, 2011(Schiel, , 2019). Nevertheless, we were able to access sites via air within a week of the earthquake to begin structured surveys of the former intertidal zone over the following few weeks when most species were readily identifiable, even down to small understory species. From these data we established the pre-earthquake condition of most of our sites, which serves as a reference for recovery over the past four years. The extensive rocky shore habitats of this coast support a rich diversity of intertidal and shallow subtidal marine species including habitat-dominating seaweeds, understory species, rock lobsters, New Zealand abalone (pâua) and other invertebrates (Lilley and Schiel, 2006;Schiel, 2006;Gerrity et al., 2020). Physical habitats include a wide variety of substrate types and topographies including near-horizontal platforms and extensive boulder fields, interspersed with dynamic mixed sand-gravel and sandy beaches.

Intertidal Community Surveys
Intertidal surveys were done at eight locations. Replicate sites (usually two) within each location were separated by at least 500 m. Sites and locations were selected to encompass the length of the earthquake-affected coastline and represent different levels of uplift. These were: control (no uplift); low uplift (0.5-1.4 m); medium uplift (1.5-2.5 m); high uplift (4-5.5 m, Figure 1). Many were areas of particular ecological, cultural, or commercial importance. Note that there were no unaffected sites to the north that could serve as "controls" due to the predominance of sandy beaches in this area. Rocky reefs and boulders occupy around a third of the earthquake-affected coastline (Gerrity et al., 2020;Orchard et al., 2021). Around 2 km of coastline experienced uplift beyond 4 m.
Surveys were done along 30 m permanent transects, one each within three intertidal zones (low, mid, and high) that are associated with characteristic flora and fauna and have been the subject of many years of surveys and experimental studies on South Island shores (Schiel, 2011). We tried to sample on the lowest tides, generally 0-0.2 above Lowest Astronomical Tide (LAT), with the tidal range being 2 m (Land Information New Zealand, 2018). Tidal heights for the three zones varied by wave exposure but were generally around 0-0.5 m, 0.6-1.2 m, and > 1.2 m above LAT. Initial sampling was completed immediately after the earthquake in the former high, mid and low zones. Most of the uplifted organisms died and disappeared within several weeks. Subsequent "post-earthquake" surveys were done in the newly formed equivalent zones, and these comprise the time series on which recovery trajectories are based. Algae and invertebrates were identified to species level when feasible or to the finest possible taxonomic resolution. Their abundances were recorded in ten 1 m 2 quadrats located randomly along the 30 m transect in each tidal zone. Abundances were expressed as% cover for sessile organisms and as counts for mobile animals (mostly gastropods).
Intertidal community structure was analyzed statistically using a distance-based permutational analysis (PERMANOVA), testing for differences on the final sampling date among uplift groups and sites within them (Uplift-fixed, 4 levels: Control, low, medium, high, and Site-random, nested within Uplift, 16 levels). Data were square-root transformed to de-emphasize the influence of abundant taxa, and analyses were based on Bray-Curtis similarities. For the Bray-Curtis similarity matrices, a dummy variable of 0.01 was used so that double zero data FIGURE 1 | The sites used for repeated monitoring were within 8 "locations," across c.130 km of coastline. The numbers in brackets indicate the number of sites per location. Uplift categories are indicated. Inset shows the km of coastline within different substratum categories (from Orchard et al., 2021). Sampling was done only in "rocky reef" and "boulder" habitats.
were treated as 100% similar. SIMPER was used to identify taxa contributing to dissimilarity between communities (PRIMER, Clarke and Gorley, 2015). The results included here mainly relate to broad taxonomic groups (i.e., groups of species sharing common morphological and life-history traits: canopyforming large brown algae, lower lying fleshy red algae, and limpets) and not to individual species. PERMANOVA was also used to test uplift and site effects on the full community in the last survey. Principal Coordinate Analysis was used for multidimensional scaling graphs.
The relationship between the degree of uplift and the% cover of large brown algae was tested using an exponential decay model with log-data. The relationship in abundance trends of large brown and fleshy red algae was further explored using a mixedeffects linear regression model including the cover of large brown algae as the response variable and that of fleshy red algae as a fixed effect. Sites were treated as random effects to account for the geographical heterogeneity of the surveys and to partition among-and within-site variability. A random intercept model was used and conditional and marginal pseudo-coefficients of determination were calculated to account for the proportion of variance explained by the fixed factor alone and the proportion of variance explained by both the fixed and random factors (so it accounts for site-by-site differences in addition to the effect of the fixed factor). To address the large variability in the data, a quantile regression model was also developed to test whether the relationship between these algal classes remained stable across the entire range of the data set. One site (Waipapa Bay 1) was excluded from analyses because large brown and fleshy red algae were absent.

Temperature and Light
HOBO temperature data loggers were placed in the low, medium and high tidal zones of each site, and maintained over the study. These were used to assess thermal conditions over the study, and were placed in areas adjacent to our survey transects. Some loggers failed (as is typical) so data sets are not always complete. Each logger was attached to a cage on the reef, so they were not in touch with the reef itself. This provided an accurate assessment of near-reef temperatures of air and water as the tide came in, but not of the surface temperatures of the rocks themselves on which organisms settled.
To gauge the light environment, PAR (photosynthetically active radiation) loggers (Odyssey R , Dataflow Systems Ltd) were deployed in the shallow subtidal zone, and above maximum high water at four locations 16 months post-earthquake (February 2018). These loggers were set to record integrated irradiance every 10 min and were downloaded and cleaned every 2 months for 1 year. Antifouling paint was applied to the top of the sensors (avoiding the sensor) and very little fouling of the sensors was observed over the course of the deployments. Using PAR data collected in the subtidal and at the surface, the daily PFD (photo flux density) was calculated. With daily PFD from the surface and at depth (average depth of the subtidal sensor accounting for tidal flux) the diffuse attenuation coefficient [Kd (PAR), hereafter referred to as Kd] was calculated daily. The daily Kd, the surface PFD and the compensating irradiance of Durvillaea antarctica/poha (Tait et al., 2015) were used to calculate the maximum habitable depth threshold for this important foundation species (see Tait, 2019). We then examined the proportion of days for which maximum habitat depth thresholds were shallower than 10 m (as an indication of light levels which may affect productivity of this species) and shallower than 5 m (as an indication of light levels which threaten growth, reproduction and survival).

RESULTS
The most obvious change to the habitats after the earthquake involved widespread algal mortality on the uplifted reefs. For example, Wairepo Reef (Kaikôura) is one of the most studied reefs in NZ and was noted for its lush beds of the desiccationresistant fucoid alga Hormosira banksii ( Table 1) across large stretches of the reef platform. Numerous experimental studies showed the high diversity of its understory, comprised of over 100 species, many of which relied on the moist canopy cover of Hormosira (e.g., Lilley and Schiel, 2006). This reef was lifted by c 0.9 m (Orchard et al., 2021). Temperatures spiked quickly after the earthquake, reaching well over 30 • C on many daytime low tides for the next month (Figure 2A). Water still covered this reef at high tide, but emersion times increased to 4-4.5 h in the semi-diurnal tidal cycle. Virtually all algae disappeared from this extensive series of platforms over several weeks . This pattern of high temperatures in the new mid and low tidal zones persisted in subsequent summers. From winter 2018 to mid-summer 2019, for example, the temperatures in these tidal zones frequently exceeded 35 • C for protracted periods from December onward ( Figure 2B). Any large algae that had recruited over the cooler months could be seen to be desiccating and then dying. This was especially evident in the summer of 2017-18 ( Figure 3A) when a severe marine heat wave affected the coast of southern NZ, and in conjunction with very low tides and hot air temperatures caused mass mortality of intertidal seaweeds along the east coast of the South Island . Experimental data showed that when a canopy of Hormosira was artificially placed onto a reef, temperatures below the canopy never reached a lethal level over the hottest part of summer , indicating the potential facilitative effects of a canopy, should one become established, on species below ( Figure 3B).
One of the formerly dominant species in the low tidal zone along the exposed coast was Durvillaea antarctica/poha (these two species occupy the same habitats and are generally indistinguishable in the field; Fraser et al., 2009Fraser et al., , 2012. Durvillaea species have a massive holdfast, thick stipe, and long leathery blades that both dampen the swell (Hay and South, 1979) and affect the understory species' composition (Santelices et al., 1980;Westermeier et al., 1994). Durvillaea had high mortality on virtually all reefs due to uplift and desiccation, and population recovery has been slow. The analysis of the maximum habitable depth threshold for Durvillaea showed that some of the most uplifted sites were frequently light-limited ( Figure 4A). For example, Waipapa Bay sites experienced > 5 m of uplift that was followed by numerous high-sediment events during rainfall because of the erosion of earthquake-damaged hills in nearby catchments, and the light environment was frequently below its compensation point (Figure 4A). An initial cover of bull kelp of around 55% at these sites has remained at zero after 4 years. There was some recovery of Durvillaea at Omihi (from around 30% cover pre-earthquake to around 15% cover after 4 years), a site of around 1.7 m of uplift. This site frequently experienced Abundant at most sites pre-EQ; virtually disappeared post-EQ and poor recovery; grows to 10 m long as adults, with great biomass (>20 kg)

Durvillaea willana Shallow subtidal
Winter >10 Greatly affected by EQ; occasionally found in lowest tidal zone; fronds grow to several m; great biomass (>20 kg per individual)

Low intertidal band
Spring-summer Prob v long-lived Very tough, dense intertwined holdfasts; fronds up to 40 cm; high initial mortality in some sites, but recovered well in sites it previously occupied

Lower-mid intertidal
Spring-summer c. 7 Common and abundant pre-EQ; high initial mortality; recovery observed in some places

Lower-mid intertidal
Spring-summer c. 7 Not as abundant as C. torulosa, tends to occur slightly lower on shore; high initial mortality, some recovery observed

Cystophora scalaris
Low intertidal Spring-summer c. 7 Least abundant of the Cystophora spp.; mostly found in tide pools and found only in a few sites in any abundance post-EQ

Mid intertidal
Year round c. 7 This species was ubiquitous in the mid-intertidal zone of rocky reefs; provides extensive cover on reef platforms, biomass up to 6 kg/m 2 wet wt. Extensive post-EQ mortality and little recovery despite being the most desiccation-resistant fucoid.

Marginarialla boryana
Subtidal Spring-summer ??? Mainly a subtidal species that occurs occasionally in the lowest tidal zone on exposed shores; high initial mortality post-EQ, but recovered slightly in some places

Lessonia variegata
Subtidal Winter ???, but probably long-lived A laminarian species; mainly subtidal that occurs occasionally in the lowest tidal zone on exposed shores; high initial mortality post-EQ, but recovered slightly in some places Only species with > 1% cover at any time are shown. Estimates of reproductive periodicity and life span are derived mostly from long-term experimental studies of the authors. ??? = unknown.
Frontiers in Ecology and Evolution | www.frontiersin.org a compromised light environment below the depth threshold. Furthermore, a moderate rainfall event in April 2018 caused an immediate and dramatic reduction in light availability for all sites monitored, the effects of which lasted for several weeks, especially at Waipapa (Figure 4B). Across the full time series, the Waipapa sites experienced very shallow maximum habitable depth thresholds 30% of the time, while the two least uplifted sites had sufficient light > 70% of the time (Figure 4C). Four years post-earthquake, algal cover varied significantly between tidal zones, across uplift levels, and at sites within uplift levels. There was very limited recovery of algae in the high and mid tidal zones of uplifted reefs, with < 10% cover of encrusting coralline and ephemeral green algae, so the focus here is on the abundance of key taxa in the low tidal zones of sites across uplift levels. In the low zone, there were significant differences between uplift levels [Pseudo-F (3 , 12) = 2.56, p < 0.01] and sites within uplift levels [Pseudo-F (12 , 144) = 9.31, p < 0.01]. The abundance of large brown algae generally decreased with increasing uplift across the control and uplift groups (Figure 5). One prominent feature was the decline in brown algal cover at 12-16 months at all sites and uplift levels, which coincided with a marine heat wave and high air temperatures combined with calm sea conditions over several days. Large brown algal cover at the Oaro control sites recovered over the next 2 years ( Figure 5A). Similarly, the large brown algae cover of low uplift sites mostly recovered to pre-earthquake levels by 4 years, although there was considerable variation among sites ( Figure 5B). Medium uplift sites remained below their pre-earthquake cover, except for one site at Omihi (Figure 5C). The high uplift sites around Waipapa, which formerly had around 70% cover of large brown algae (mostly Durvillaea spp.) had < 20% cover after 4 years, with one site having no brown algal cover at all (Figure 5D). The large brown algae species found on uplifted recovering reefs were primarily Carpophyllum maschalocarpum and Cystophora spp., with some Marginariella boryana and a small cover of the kelp Lessonia variegata at a couple of sites (cf . Table 1). Overall, the low-uplift and control sites had the highest average large brown algal cover (c. 60%), followed by the medium-uplift (38%), and the high-uplift sites (8%).
Fleshy red algae were the other dominant group of habitatdefining seaweeds along the coast, and they were generally not as affected by the heat wave in the summer of 2017-18 as were large brown algae (Figure 6). After 4 years, there were significant differences between uplift levels [Pseudo-F (3 , 12) = 1.82, p < 0.05] and sites within uplift levels [Pseudo-F (12 , 144) = 11.10, p < 0.01]. At the control sites, red algae fluctuated over the years and had a slightly lower cover (-10%, Table 2) after 4 years than at the start of the study (Figure 6A). Red algae in low uplift sites generally had the same cover throughout the study and seemed to be little affected by the earthquake (Figure 6B and Table 2). However, medium uplift sites showed increases of 25% in cover (across all sites, Table 2) over the 4 years ( Figure 6C), and high uplift sites showed overall increases of 22% in red algal cover over the years in comparison to pre-earthquake conditions ( Figure 6D). The high-uplift areas had the greatest variation among sites. The taxa of fleshy red algae primarily responsible for these patterns were Gelidium microphyllum, Pterocladia lucida, Ceramium spp., Gigartina chapmanii, Chondria macrocarpa, Polysiphonia spp., Echinothamnion spp., and Champia sp.
Two prominent features of the recovery process were the relationships between uplift and percentage cover of large brown algae, and between the percentage cover of brown and fleshy red algae. There was an exponential fall-off in brown algal cover with degree of uplift (Y = 75.5 e −0.499X ; r 2 = 0.89, Figure 7A) after 4 years. Regression analyses highlighted a negative relationship between the cover of large brown algae and fleshy red algae. This relationship accounted for 10% of the variability in the abundance of large brown algae in the mixed-effects linear regression, whereas the whole model (accounting also for site-by-site differences) accounted for 36% of the variability. Quantile regression analyses confirmed the presence of a strong negative relationship between the two groups across five different quantiles (q10 slope -0.25, p < 0.001; q25 slope -0.35, P < 0.001; Q50 slope -0.41, p < 0.001; q75 slope -0.44, p < 0.001; q90 slope -0.4, p < 0.01, Figure 7B). These relationships contribute to the observed percentage gains and losses of red and brown algae across tidal zones in the different uplift levels ( Table 2).
Limpets are the dominant grazers at all study sites and were typically abundant from the mid tidal zone upwards. At 6 months post-quake, control and low uplift sites had the greatest densities of limpets (c. 30 m −2 ) but after 4 years had the fewest (< 20 m −2 ; Figure 8A). The medium and high uplift sites showed increases in limpets over the 4 years. In particular, the medium uplift sites  had an average of 60 limpets m −2 , but this was driven mainly by one of the sites (Okiwi Bay) within this uplift level that had 216 limpets m −2 . There was therefore a significant difference in sites within uplift levels [Pseudo F (12 , 144) = 2.99, p < 0.01], but not among uplift levels [Pseudo F (3 , 12) = 2.66, p = 0.10].
The limpets were mostly Cellana species, but also included Notoacmea and pulmonate limpets (Siphonaria). The decline in limpet numbers in the control and low-uplift sites did not reflect a widespread decline in the abundance of all limpet species, but was due to the absence of large clusters of Siphonaria spp. Another gastropod, the commercially valuable pâua (abalone, Haliotis iris) recruited well in all post-earthquake years, with numbers accumulating through to year 4 ( Figure 8B). This was related to recovering juvenile habitat of small boulder-fields in the lowest tidal zone, and accumulating numbers of reproductive adults in most shallow rocky habitats (Gerrity et al., 2020). Taken together, there remained large differences in community structure between sites and uplift levels within tidal zones at 4 years post-quake (Figure 9). Multivariate analysis showed uplift had a significant effect on intertidal community composition in the post-earthquake high tidal zone [Uplift: Pseudo-F (3 , 12) = 2.42, p < 0.05] as did sites within uplift levels [Pseudo-F (12 , 144) = 4.96, p < 0.01; Figure 4A]. In particular, the control and the low-uplift groups differed from the mediumuplift group because of higher covers of ephemeral green (Ulva spp.) and red algae (Pyropia spp.). However, this was a reflection of occasional blooms in the abundance of ephemeral algae rather than a real earthquake legacy. Long-lasting earthquake FIGURE 7 | (A) The relationship between percentage cover (± SE) of large brown algae across degrees of uplift. (B) Relationship between the abundance of large brown and fleshy red algae 3 years after the earthquake estimated through a quantile regression model. Regression lines are displayed for significant relationships. Symbol colors indicate different levels of uplift (white = no uplift, green = low uplift, yellow = medium uplift, red = high uplift). Data from the high-uplift site of Waipapa Bay 1, where both large brown and fleshy red algae were not present, were not included in this analysis. effects were especially evident in mid zone, where abundant algal communities occurred only at control sites ( Figure 9B). This led to a significant effect of uplift [Pseudo-F (2 , 8) = 2.62, p < 0.01], but there was also a significant site effect [Pseudo-F (8 , 99) = 9.37, p < 0.01]. In the low zone, the composition of benthic communities was different in the low-uplift group (where large brown algae, particularly Carpophyllum maschalocarpum, had recovered) compared to the medium-and high-uplift groups [which were characterized by high percentage cover of red algae; Uplift: Pseudo-F (3 , 12) = 2.43, P < 0.01, Figure 9C]. No other uplift group differed from the others but there was significant variability in the structure of benthic communities among sites within each uplift group [Pseudo-F (12 , 144) = 7.66, P < 0.01, Figure 9C].

DISCUSSION
This study has been instructive in understanding the dynamics of recovery from a major disturbance event involving species, communities and the infrastructure on which they build and depend. It was unusual, if not completely novel, to witness such extensive damage to a coastal ecosystem. Rocky reef communities that had been dynamic but intact in structure and function over many decades in the face of periodic storm and wave impacts (Schiel, 2011(Schiel, , 2019 were obliterated by a major event that occurred over a matter of minutes . The horizontal distance between the former and post-earthquake high tide marks reached c. 200 m in areas of major uplift. It was somewhat disconcerting to see entire shallow coastal assemblages high and dry out of tidal influence. This, however, provided the FIGURE 9 | Principal coordinates analysis (PCO) plots showing differences in the composition of benthic communities in the post-earthquake high (A), mid (B) and low zone (C) across sites with different degrees of uplift 48 months after the earthquake. The symbols represent the centroid of each site and the colors the different levels of uplift (white = no uplift, green = low uplift, yellow = medium uplift, red = high uplift). Sites are ordered north to south within each uplift group. Only the high and the low zone were sampled at high-uplift sites.
opportunity to assess pre-earthquake communities quantitatively that were previously difficult to access and therefore rarely studied in detail.
At most sites, uplifted reef platforms were not compensated by new rocky substrates uplifted from the subtidal zone; this led to a new configuration of low and mid tidal zones that are now near-vertical in places. This morphological change was important in recovery because there was a much smaller tidally influenced zone in which communities could re-assemble. A surprise that unfolded within a month post-earthquake was the loss of speciose algal communities, which had been the focus of decades of studies, in the mid and upper tidal zones of extensive reef platforms, even though they still had tidal coverage. Furthermore, recovery of the low tidal zone was substantially set back by the severe heat wave a year after the earthquake, the timing of which coincided with the replacement of large brown algae by a suite of tough, fleshy red algae. The density and longevity of the latter have precluded settlement of large brown algae in the period since the heat wave. Intertidal platforms and boulders, composed of Paleocene limestone overlain by Oligocene gray mudstones (Kirk, 1977), have also eroded at annual rates > 35 mm at several sites post-quake  compared to a long-term average of 1.13 mm yr −1 (Kirk, 1977;Stephenson and Kirk, 1998).
Studies from Chile have shown that recovery from coastal uplift can take several years. For example, Castilla (1988) found that uplift of 0.4-0.6 m caused extensive mortality of the laminarian alga Lessonia nigrescens, which had formed a conspicuous band on the low shore. As understory species also died, the areas they occupied became covered in barnacles, and Lessonia populations were still recovering 3 years after the earthquake. In 2010, a mega-earthquake of Mw 8.8 caused coastal lifting up to 3.1 m (Castilla et al., 2010). Similar to our study, there was widespread mortality of intertidal organisms after being exposed to solar radiation and elevated temperatures, especially Durvillaea antarctica and red algae which became bleached and desiccated. Jaramillo et al. (2012) also documented broad-scale impacts of this earthquake and found, as in our study, that the ecological impacts varied strongly with magnitude and direction of land-level change across different shore types and with the mobility of characteristic biota. Because of structural changes to the intertidal zones, the likelihood of recovery to a former state was unclear.
One of the novel features of recovery of the shores in southern NZ is the role of mostly long-lived fucoid algae that are the dominant habitat-formers. These generally have short dispersal distances of propagules, often a matter of meters, and rely on drifting, reproductively active adults for long distance dispersal (Schiel, 2011). Because of the patchy nature of post-earthquake fucoid populations, connectivity to previously occupied areas was compromised and mostly relied on drifting adults for recolonization. This is illustrated by one of the hardest hit taxa, Durvillaea. Chilean studies have shown that drift Durvillaea antarctica can be abundant onshore (Tala et al., 2019) and it is known that this species can drift extensively in southern seas (Waters, 2008). However, the deposition of drifting fucoids at distant inshore sites is probabilistically low over short time periods (Hawes et al., 2017). Additionally, Durvillaea is dioecious, so male and female fronds must arrive inshore together, and this needs to coincide with their relatively short reproductive season (c. 8 weeks) in winter. Few drift Durvillaea have been seen at our study sites since the earthquake. As well, its sporelings are known to be highly vulnerable to heat stress (Hay, 1979), and in at least some sites of former abundance the water clarity is so poor that effective growth would be compromised.
Similar impediments to establishment are faced by other fucoids. For example, Hormosira banksii is the most desiccationresistant fucoid because of its mucilage-filled fronds (Brown, 1987). However, the reduced period of tidal immersion, erosion of reefs and increases in fine sediments can all compromise effective recruitment (Alestra and Schiel, 2015). This is coupled with low-tide temperatures frequently exceeding 35 • C during summer, killing off the sparse annual recruitment in the lower mid-tide zone and effectively turning this perennial and formerly dominating species into an ephemeral one. This species has declining productivity beyond 25 • C (Tait and Schiel, 2013), and its canopy interacts synergistically with other fucoids, such as Cystophora torulosa and Carpophyllum maschalocarpum, to increase per-area primary productivity (Tait andSchiel, 2011, 2013), which is similar to fucoids elsewhere (Colvard et al., 2014). By comparison to the low-shore fucoids, fleshy red algae seem to be more resistant to heat stress and can remain productive and recover from desiccation at temperatures above 30 • C (Smith and Berry, 1986).
In contrast to algae, broadcast-spawning invertebrates recovered quickly, as would be expected. In the first 6 months post-earthquake, the mid and lower intertidal zone was bright green along much of the coast from a massive bloom of mostly Ulva spp. As grazers recruited, such blooms became more patchy and ephemeral. One surprise, however, was the recovery of pâua (abalone, Haliotis iris) populations in the low intertidal-subtidal margins. Despite the initial loss of recruitment habitat in the lower zone, there was good recruitment in each year following the earthquake (Gerrity et al., 2020). Adult pâua were lifted into this zone from uplifted subtidal reefs and also migrated there from deeper areas, and there was a ban imposed on recreational and commercial fishing along the coast. The combined result was an abundance of shallow pâua populations at levels not seen in decades.
Resilience and recovery clearly need to be considered within the context of species' life histories, life spans, their interactions with each other, and prevailing conditions, which can involve long time trajectories. In this study, recovery has been influenced by a combination of altered topography, increased temperatures at the reef surface, erosion of rocky substrata, sedimentation from earthquake-damaged catchments, a compromised light environment in nearshore waters, and human-induced stressors from increased coastal access around headlands. These processes are acting on a pattern of widespread initial mortality that occurred even at lesser uplifted sites, and are consistent with major re-assembly processes being driven by relatively small increments of relative sea-level change (Orchard et al., 2021). The variable patterns of recovery we have observed reflect the influence of complex site-specific combinations of stressors that are mostly hindering the re-establishment of supportive ecological infrastructure. Exceptions, however, include the recovery of large brown algae at several of the low uplift (<1 m) sites which likely reflects the influence of surviving population remnants, and the expansion of fleshy red algal cover at others in absence of the recovery of brown algae. In relation to the degree of uplift, there are confounding factors among sites from the rate of weathering characteristic of softer substrata that predominate at them. The associated erosion effects are a likely cause of slow recovery on these reefs and exhibit further interactions with turbidity and alterations to the light environment that we recorded in these areas.
It may turn out that there are multiple stable states depending on the level of uplift and its interactions with other factors. These complex interactions affect how communities re-align themselves over longer time periods but notably exert their most pervasive effects through the key habitat-forming species that generate the ecological infrastructure and supportive conditions upon which others depend. Scheffer et al. (2001) point out that state shifts in ecosystems can be reflected more by biological than physical factors because biotic feedbacks can stabilize communities in different ways, and that such shifts may be triggered by catastrophic physical events. So far, our time scales for recovery have not coincided with those needed for stabilization of the most earthquake-affected ecological infrastructure, notably that provided by the large brown fucoid species characteristic of this coast.
Much of the resilience and ecological infrastructure literature addresses nature's goods and services, and potential management interventions to maintain them (da Silva and Wheeler, 2017; Chambers et al., 2019). It would initially seem that there would be few management strategies for ecosystem recovery that could be effective in the face of cataclysmic earthquake impacts followed by a major heatwave. And yet manageable stressors and potentially effective pathways to facilitate recovery trajectories can be identified. For example, there is a current legislative initiative to restrict vehicular traffic on the recovering beaches and tidal platforms around Cape Campbell in the north of the study area. The coast-wide closure of the pâua fishery also resulted in an unanticipated and quick rebound of pâua stocks, which led to current proposals favoring the reduction of commercial and recreational fisheries, and an adoption of a more adaptive and resilience-based management strategy. Additionally, there were lightly affected areas representing ecological "safe havens" along the coast at locations such as Oaro, where traditional management practices of restricted fishing have been overseen by customary Mâori tikanga (practices). These highlight the value of multiple protected areas along coastlines which, perhaps by chance and good fortune, can act as insurance policies to provide sources of reproductive propagules to aid the recovery of key habitat-forming species. These examples illustrate that a focus on the protection of ecological infrastructure is a tractable objective for the facilitation of disaster recovery. It can be promoted by the strategy of addressing manageable stressors such as sediment loss and overharvesting, and is indeed essential for the return of desirable states in the aftermath of natural disasters.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
DS wrote the manuscript and led the research program. SG and TA led the field team and collected much of the data. TA and TF did data analyses. RD led subtidal research and helped set up the program and did initial surveys. SO led community liaison and did the analyses of physical changes to the coastline. MT analyzed heat wave effects. LT did photophysiology and light studies. All authors contributed to the article and approved the submitted version.

FUNDING
This research was funded by the MBIE contract UOCX1704 and the Ministry of Primary Industries grant KAI2016-05. Aligned funding was by the Sustainable Seas National Science Challenge (CO1X1901).

ACKNOWLEDGMENTS
We thank Dr. John Pirker for iwi liaison, Jason Ruawai and Te Runanga o Kaikoura for input and support throughout the study, Dr. Sharyn Goldstien for helping with logistics and involvement in the hectic days following the earthquake, NCTIR for granting us access to field sites during the extensive period or cliff stabilization and road closure, Dr. Rich Ford of MPI for support throughout the program, Dr. Mike Hickford for help with logistics, Dan Crossett who helped with field work in parts of the program not reported here, and many others who helped at times with field work and other support.