The ocean is experiencing long-term warming due to anthropogenic climate change, having absorbed over 90% of the excess heat energy produced by increased greenhouse gases (IPCC, 2021). This long-term increase in ocean temperature is the main factor driving an increase in the frequency, extremity and duration of marine heatwaves (MHWs; Fragkopoulou et al., 2023; Schlegel et al., 2019; Smith et al., 2021), defined as more than five days of ocean temperatures above the 90^th percentile of historical observations (Hobday et al., 2016). MHWs have immediate and significant impacts on marine ecosystems (Frölicher et al., 2018; Smith et al., 2021). Consequences for fisheries include biomass reductions (Cavole et al., 2016; Cheung and Frölicher, 2020; Gammelsrød et al., 1998; Mills et al., 2013; Ñiquen and Bouchon, 2004), which can be much greater and more rapid than those caused by decadal-scale warming (Baudron et al., 2014; Cheung and Frölicher, 2020; Givan et al., 2018).

In Aotearoa New Zealand, wild-capture fisheries exports were worth over NZ$1.5 billion in 2023 (MPI, 2023), with inshore finfish contributing approximately $75 million (MPI, 2023). In addition, many coastal species also hold important recreational and cultural value (Bess, 2016; McCarthy et al., 2013, 2014; Parsons et al., 2014). Future increases in ocean temperature are predicted to reduce growth rates of key New Zealand species, such as snapper, Chrysophrys auratus, and tarakihi, Nemadactylus macropterus (Cummings et al., 2021; Morrongiello et al., 2021), and the recent 2017–2018 summer MHW caused mass mortalities of caged coastal aquaculture salmon (Salinger et al., 2019). However, wild fishes can track their preferred thermal niche, unlike sessile species that can only do so on a population-turnover time scale, or aquaculture fish whose range is limited to the cage frame. Many fish species have undergone poleward range shifts in response to warming (García Molinos et al., 2016; Last et al., 2011; Parmesan and Yohe, 2003; Pinsky et al., 2013; Poloczanska et al., 2013), including in New Zealand (Lavin et al., 2023; Middleton, 2022).

Range shifts are not limited to latitudinal movement, however, as the ocean is a three-dimensional space and temperatures decrease with depth (Bongaerts et al., 2010; Frade et al., 2018; Leichter et al., 2006; Riegl and Piller, 2003). Several fish species have demonstrated long-term deepening of their distributions in response to warming (Chaikin et al., 2021; Dulvy et al., 2008; Givan et al., 2018; Nye et al., 2009; Perry et al., 2005; Rutterford et al., 2015), and deeper waters may also offer short-term refuge from thermal stress. Mesophotic ecosystems extend from where light levels are approximately 1% of surface irradiance, up to a maximum depth beyond which primary production cannot occur (Cerrano et al., 2019). These mesophotic conditions can occur across a broad range of depths depending on local conditions but are often found between ~30 and 150 m (Cerrano et al., 2019). The deep reef refuge hypothesis (DRRH), proposed originally for mesophotic coral ecosystems, suggests that mesophotic reefs may be buffered against disturbances that impact the shallows, and may be able to re-seed shallow communities following disturbance (Bongaerts et al., 2010). It is now understood that the DRRH holds only for some coral species under certain disturbances (Bongaerts et al., 2017), but its role in supporting mobile organisms such as fish is less well-studied. Mesophotic reefs could act as ecologically viable thermal refuges from MHWs, particularly for coastal reef fishes that rely on benthic habitat for structural complexity and food. One of the few studies exploring this showed that the subtropical reef fish Lethrinus miniatus on Australia’s Great Barrier Reef demonstrated short-term movement from euphotic to mesophotic reefs when mean daily water temperatures increased (Currey et al., 2015). The potential for a deep reef refuge has not yet been explored for temperate coastal fishes.

While the majority of recent physical oceanography studies found the intensity and duration of MHWs increase with depth (Dayan et al., 2023; Fragkopoulou et al., 2023; Hu et al., 2021; Schaeffer et al., 2023; Schaeffer and Roughan, 2017; Sun et al., 2023; Zhang et al., 2023), several studies investigating the ecological impacts of MHWs across depth have indicated a mesophotic thermal refuge effect, including for canopy-forming macroalgae (Giraldo-Ospina et al., 2020), corals (Frade et al., 2018; Muir et al., 2017), sponges (Perkins et al., 2022) and anemones (Haguenauer et al., 2021). Physical oceanography and ecological impact research into mesophotic MHWs are still somewhat divorced; physical oceanography studies have not yet analysed actual or modelled MHW impacts across depths, and though some ecological studies include temperature-depth data (Frade et al., 2018; Haguenauer et al., 2021), none have yet conducted MHW analyses on those data. Furthermore, to our knowledge, all subsurface MHW studies to date focus on heatwaves according to depth-specific temperature regimes. This is valuable when considering impacts on species which reside within specific depth strata and are acclimatised to their usual temperature regimes, particularly sessile species. However, for mobile organisms with preferred depth ranges, the buffering potential of deeper depths against MHWs (relative to their preferred depth range’s temperature regime) is highly relevant. In addition, depth-driven differences in the maximum absolute temperatures that organisms are exposed to during MHWs might be a more important predictor of stress and physiological response (Cavanaugh et al., 2019).

Exploring the ecosystem-level implications of mesophotic reefs as thermal refuges from MHWs requires consideration of ecosystem dynamics and species interactions. Multispecies size-spectrum models that incorporate temperature effects on physiological rates are well suited for this task. Multispecies size-spectrum models (MSSMs) incorporate species-specific life history traits and well-evidenced size-based feeding theory (Blanchard et al., 2009). They allow for ontogenetic diet shifts, opportunistic predation and cannibalism (Jennings et al., 2001), whilst requiring relatively few parameters compared to data-hungry end-to-end models such as Atlantis (Fulton et al., 2004, 2005, 2011) and Ecopath with Ecosim (Christensen and Walters, 2004). In recent years, especially following the development of the R package mizer (Delius et al., 2023; Scott et al., 2014), MSSMs have been used to investigate fisheries management strategies (Benoit et al., 2022; Blanchard et al., 2014; Robinson et al., 2022; Spence et al., 2021; Szuwalski et al., 2017), blue carbon stocks (Falciani et al., 2022), stock-recruitment relationships (Canales et al., 2020), ecosystem recovery signals and ecological indicators (Clements et al., 2019; Zhang et al., 2018), fisheries impacts (Forestier et al., 2020), and the effects of seasonal processes on dynamics and yields (Datta and Blanchard, 2016). Woodworth-Jefcoats et al. (2019) simulated the impacts of ocean warming and fishing on the bigeye tuna longline fishery in Hawai’i and found that detrimental warming effects could be mediated but not completely mitigated by fisheries management changes. These findings were supported by an Ecopath with Ecosim model of the same system (Howell et al., 2013). Audzijonyte et al. (2022) included size-structured benthic resources as well as temperature effects for a Tasmanian coastal rocky reef ecosystem and found that temperature-driven changes in planktonic and benthic resources affected fish biomasses and yields more strongly than temperature-driven physiological changes did. Lindmark et al. (2022) modelled temperature effects on planktonic resource spectrum dynamics as well as on fish physiology and found that warming resulted in increases in size-at-age and growth rates but decreases in large fish abundance and fisheries yields. In investigating similar warming impacts, Kuo et al. (2022) also incorporated species distribution models, allowing species to respond spatially with horizontal range shifts, habitat shrinkage and changes in species overlap and interactions.

Despite this growing body of research on temperature effects in ecosystem models, the potential of deeper waters to act as a thermal refuge has not yet been explored, although it has been highlighted as an important research avenue (Audzijonyte et al., 2022; Woodworth-Jefcoats et al., 2019). In this study, we use modelled temperature-depth data and develop a temperature-dependent MSSM to test the hypotheses that in temperate ecosystems, (1) MHWs are buffered at mesophotic depths compared to euphotic depths, and (2) that mesophotic thermal refuge use mitigates reductions in fish biomass and productivity during MHWs in temperate coastal ecosystems.

Mesophotic reefs at the Poor Knights Islands have the potential to act as refuges from MHWs for mobile species. During summer, when temperatures are most likely to approach the upper bounds of coastal fish species’ thermal tolerance ranges (Shultz et al., 2016), mesophotic depths were cooler than euphotic depths; for example, 50 m was 3.22°C cooler than 10 m. Temperatures were more vertically mixed during winter, consistent with the typical seasonal pattern observed in shelf seas (Frade et al., 2018; Rippeth, 2005; Sallée et al., 2021; Simpson and Rippeth, 1993). Importantly, the duration, intensity and frequency of MHWs, which are considered the most dangerous present and future source of heat stress for marine species (Bell et al., 2023; Frölicher et al., 2018), were buffered from 30 m depth and below at the Poor Knights relative to MHWs at euphotic depths ( Figure 4 ). In contrast, recent studies using depth-specific MHW thresholds have found that subsurface MHWs are only occasionally less intense at mesophotic depths (Zhang et al., 2023) and are often equally or more intense and have longer durations than at euphotic depths (Amaya et al., 2023; Elzahaby and Schaeffer, 2019; Fragkopoulou et al., 2023; Großelindemann et al., 2022; Hu et al., 2021; Scannell et al., 2020; Schaeffer et al., 2023; Schaeffer and Roughan, 2017; Zhang et al., 2023). Part of the explanation for these contrasting results is due to our definition of the MHW, which is determined by threshold temperatures at euphotic, shallow depths. However, differences in scale and the oceanographic processes at play compared to other studies may also play a role in this result. Mesophotic reefs at the Poor Knights potentially receive a regular (twice-daily) input of cool water at depth due to storms at sea and shoaling internal waves and/or bores occurring at tidal frequencies, which are a well-known feature of the northeast continental shelf of New Zealand (Sharples et al., 2001; Sharples and Zeldis, 2019) and the Poor Knights Islands specifically (Stevens et al., 2005). The Moana Hindcast, like similar regional ocean configurations using ROMS, is expected to resolve the cooling associated with internal waves (Suanda et al., 2017; Fagundes et al., 2020). This process is unlikely to be relevant to studies such as Elzahaby and Schaeffer (2019); Fragkopoulou et al. (2023); Großelindemann et al. (2022); Hu et al. (2021); Scannell et al. (2020), or Zhang et al. (2023), all of which focus on MHWs in the open, deep (>1000 m) ocean. Shoaling internal waves are, however, well-recognised as a process that buffers heat exposure at depths similar to those considered here on coral reefs (Storlazzi et al., 2020; Wyatt et al., 2020), and may be a key process explaining the buffering effect of MHWs at mesophotic depths at the Poor Knights Islands. Further investigation of this possibility would require dedicated observations and/or numerical modelling that are beyond the scope of the present research.

Some studies of the ecological impacts of MHWs by depth – as opposed to solely oceanographic studies – have found that sessile mesophotic species, having acclimatised to mesophotic temperature regimes, experience similar thermal stress to their euphotic counterparts (Bongaerts and Smith, 2019; Smith et al., 2016: Bell et al., 2024). However, often MHW impacts on sessile species at mesophotic depths are less than those at euphotic depths (Frade et al., 2018; Giraldo-Ospina et al., 2020; Haguenauer et al., 2021; Muir et al., 2017; Perkins et al., 2022; Bell et al., 2024). This, combined with our findings of positive mesophotic thermal refuge effects for mobile fishes, warrants the consideration of incorporating assessments of some measure of euphotic- or surface-relative temperature anomalies or absolute temperatures alongside depth-specific temperature anomalies in subsurface MHW analyses that intend to predict ecological impacts.

Under temperature-dependent multispecies size-spectrum model simulations of various heatwave scenarios, most species and guilds responded with progressively larger biomass reductions with each increase in heatwave category level, but the thermal refuge effect mitigated these reductions. For fisheries species, average biomass losses were reversed (28% difference) during the moderate heatwave, negated (24% difference) during the strong heatwave, and mitigated by 21% and 20% during the severe and extreme heatwaves. The responses of herbivores, large piscivores and cryptobenthic fishes varied, however, with mesophotic refuges having a positive, neutral, or negative impact on biomass depending on heatwave category level; when the effects were positive it was likely because the temperatures of a given heatwave resulted in fewer predators and/or competitors, while not yet exceeding their own thermal tolerance limits. Population productivity, a rate which is more relevant to fisheries than biomass alone, increased slightly or was similar during historic daily average euphotic temperatures for all species during the moderate euphotic heatwave, likely due to faster growth (Carozza et al., 2018; Pepin, 1991). Temperature increases have been shown to result in there being smaller, faster-growing fishes in a population (Baudron et al., 2014). However, productivity decreased during strong, severe, and extreme heatwaves. For fisheries species, moving to a mesophotic thermal refuge almost negated average productivity losses during the strong heatwave (5% difference), negated them during the severe heatwave (17% difference) and mitigated them by 19% during the extreme heatwave.

These findings emphasise the well-evidenced negative impacts of thermal stress on fish population biomass and productivity (Capitani et al., 2022; Carozza et al., 2018; Free et al., 2019; Wilson et al., 2021), as well as the importance of species interactions in determining responses to warming (Baum and Worm, 2009; Kirby and Beaugrand, 2009; Steneck, 2012) and thermal refuge availability. Reductions in biomass and fisheries yields with ocean warming have been observed using other temperature-dependent mizer models (Audzijonyte et al., 2022; Kuo et al., 2022; Lindmark et al., 2022; Woodworth-Jefcoats et al., 2019), and Kuo et al. (2022) highlighted evidence of warming-induced trophic cascades as changes in spatial overlap and decreases in fish biomass resulted in species having to feed more on background resources. However, the community-wide negative response to the extreme euphotic heatwave highlights that the impacts of approaching or exceeding thermal tolerance limits outweigh the benefits of release from competition and predation. Of particular interest is the fact that a mesophotic thermal refuge mitigates even the effects of this year-long extreme heatwave for all species, given that by 2100 under RCP8.5, marine ecosystems will be under constant heatwave conditions relative to current thresholds and the majority of events will be extreme (Oliver et al., 2019).

This study demonstrates the utility and potential of using multispecies size-spectrum models to explore the impacts of thermal refuges on fish community biomass and productivity under climate change scenarios, but there are several features that would increase the accuracy of this size-spectrum model’s representation of a coastal rocky reef ecosystem that could be incorporated in future research. Mizer models in general are currently particularly well-suited to modelling open-ocean pelagic ecosystems, but developments have been occurring for realistically representing coastal ecosystems. Rather than having herbivores feed on the plankton spectrum, algae could be modelled explicitly and adjustments could be made to therMizer’s mechanics to incorporate the feeding on benthic resources allowed for in the mizer extension package mizerMR (Audzijonyte et al., 2022). The response of herbivores in the simulations run using this study’s model should be treated with some caution. Though the effects of temperature and predation were represented in herbivores’ response, so too was the effect of competition for the plankton resource, which is not applicable to them in a real system. Temperature effects on the planktonic (Woodworth-Jefcoats et al., 2019) and benthic resource spectra (Audzijonyte et al., 2022) could also be included, which would allow for a more complete picture of heatwave impacts. How to realistically incorporate large transient predators such as sharks, which do not act within model boundaries and whose life history traits are very different to those of teleost fishes (and hence were excluded from the model presented here), is also a valuable future research avenue.

While the observed abundances to which this size-spectrum model was calibrated were from euphotic depths, the depth ranges of a number of species captured in the model do already include upper mesophotic depths (Allen, 1991; Heemstra and Randall, 1993; Fisheries New Zealand, 2023; McMillan et al., 2011; Mundy, 2005; Paulin, 1990; Randall, 1999; Stevenson, 2004). Future models could incorporate existing depth ranges, and then allow individuals in the upper portions of the depth ranges to deepen during MHWs while individuals already in the mesophotic would remain there. However, this would require fish survey data at mesophotic depths, which are beyond recreational diving limits and require either technical divers, remotely operated vehicles or autonomous underwater vehicles to survey (Armstrong et al., 2019; Munoz et al., 2017; Pyle, 2019; Turner et al., 2018). In mizer models of fished systems, biomasses at depth are usually taken from catch data (Blanchard et al., 2014; Jacobsen et al., 2017; Szuwalski et al., 2017; Woodworth-Jefcoats et al., 2019), and it would be interesting to explore the impact of mesophotic thermal refuge availability on fisheries yields specifically. In models encompassing large areas, it would also be pertinent to incorporate latitudinal range shifts as a simulation possibility (Kuo et al., 2022). Furthermore, site-attached species such as triplefins and eels (Abrams et al., 1983; Nagelkerken et al., 2016) would be less likely to move to a mesophotic reef during warming, even if one were available. Herbivores would likely only be able to utilise the upper limit of the mesophotic as a thermal refuge, as their food is light limited; macroalgal cover reduces from approximately 15% at 30 m to 2% at 40 m and below at the Poor Knights (Harris et al., 2021). Other factors that could limit the vertical migration of different fish species include the effects of density and prior residency, and potential local niche saturation (Pinheiro et al., 2023). With some species remaining and some escaping, species interactions would likely play an even larger role in thermal refuge availability responses than they did here. However, as interesting, complex and more accurate mizer models incorporating transient predators, algal and benthic resources, existing depth ranges, latitudinal range shifts, complex species interactions, and temperature impacts on planktonic, algal and benthic resources would be, such models are beyond the scope of this study. There is also the principle of parsimony and the issue of requiring sufficient data to accurately model all these processes and interactions. Simpler models do have the benefit of being able to clearly and succinctly explore one main research question, especially when data are available but not extensive.

This study provides clear evidence supporting the hypothesis that mesophotic reefs can act as an effective thermal refuge from MHWs for commercially and recreationally important coastal New Zealand fish species. Species that responded positively to heatwave events in this study were almost always non-fisheries species, while those that suffered the greatest losses were almost always species targeted by fisheries. If these fisheries species do in fact use mesophotic refuges during marine heatwaves, it would heighten their resilience, and support the sustainability of fisheries in the face of climate change. This has significant conservation and management implications, as mesophotic reefs are often only incidentally protected when they happen to occur in marine protected areas (Bell et al., 2022). Though often buffered from thermal stress, as shown here and by several other studies (Giraldo-Ospina et al., 2020; Haguenauer et al., 2021; Muir et al., 2017; Perkins et al., 2022), mesophotic reefs are not invulnerable to anthropogenic impacts (Bell et al., 2022; Frade et al., 2018; Rocha et al., 2018), and targeting mesophotic reefs specifically for protection could improve the future climate resilience of New Zealand’s coastal fisheries.