Loss of Giant Kelp, Macrocystis pyrifera, Driven by Marine Heatwaves and Exacerbated by Poor Water Clarity in New Zealand

Marine heatwaves (MHW) are becoming stronger and more frequent across the globe. MHWs affect the thermal physiology of all biological organisms, but wider ecosystem effects are particularly impactful when large habitat-forming foundation species such as kelps are affected. Many studies on impacts from MHWs on kelps have focused on temperature effects in isolation, except for a few studies that have integrated co-occurring stress from grazers, wave exposure and nutrient limitation. It is likely that many stressors act in concert with MHWs and exacerbate their effects. Here we analyzed satellite images over 60 months to assess temporal changes in abundance of surface canopies of the giant kelp Macrocystis pyrifera in the New Zealand coastal zone. The analysis encompassed the most extreme MHW on record (2017/18), across a 6° latitudinal gradient of four regions southward from the northern distributional limit of Macrocystis along mainland New Zealand. We tested the association of surface canopy cover of Macrocystis with sea surface temperature, temperature anomalies, chlorophyll-a (a proxy for nutrient availability) and water clarity (diffuse attenuation coefficient). We found a reduced cover of Macrocystis across all regions during and after the 2017/18 MHW, with least impact at the most southern region where the maximum temperatures did not exceed 18°C. There was also an important and significant interaction between temperature and water clarity, showing that temperature-induced kelp loss was greater when water clarity was poor. These results show that notable negative effects occurred across the coastal range of this foundation species and highlight the importance of studying MHW effects across latitudinal gradients and in concert with other co-occurring stressors.


INTRODUCTION
Warming ocean temperatures, temperature anomalies, and marine heatwaves (MHW) are becoming stronger, longer and more frequent (Frölicher et al., 2018;Hobday et al., 2018;Gupta et al., 2020). Temperature affects all aspects of biology, from subcellular biochemical reaction rates to reproductive success, evolutionary mutation rates and selection pressures (Hoegh-Guldberg and Bruno, 2010;Doney et al., 2011). It is therefore not surprising that temperature anomalies can have pervasive impacts on marine productivity (Beaugrand and Reid, 2003;Montie et al., 2020), community structure McPherson et al., 2021) and ecosystem processes (Smale et al., 2019;Gupta et al., 2020). These ecological changes are particularly noticeable when warm anomalies and MHWs affect large habitat-forming foundations species, like seagrass (Thomson et al., 2015;Arias-Ortiz et al., 2018), coral reefs (Le Nohaïc et al., 2017) and kelp (Wernberg et al., , 2016. Impacts from MHWs on kelp have been documented in Tasmania , Western Australia (Wernberg et al., , 2016, North America (Arafeh-Dalmau et al., 2019;Rogers-Bennett and Catton, 2019) and Europe , and have, in extreme cases, caused range contractions and regional extinctions Wernberg et al., 2016;Straub et al., 2019;Smale, 2020). High temperatures often have predictable physiological effects on kelp tissue and population growth rates (Tait, 2014;Mabin et al., 2019). However, impacts on abundances, communities and ecosystem functioning are complicated by indirect effects, altered competitive hierarchies, the timing and magnitude of specific MHWs as well as latitudinal distribution and thermal responses of individual kelp species (Doney et al., 2011;Wernberg et al., 2016;Smale et al., 2019). For example, it has been shown that MHWs can have a strong negative impact on kelps near their warm equatorial range limit, but no or positive effects near their colder poleward range (Wernberg et al., 2010(Wernberg et al., , 2016Reed et al., 2016;Arafeh-Dalmau et al., 2019;Cavanaugh et al., 2019), highlighting the importance of measuring impacts and variation across latitudes.
To date, most research about impacts from MHW on kelp forests have focused mainly on temperature effects alone (Reed et al., 2016;Wernberg et al., 2016;Arafeh-Dalmau et al., 2019;Cavanaugh et al., 2019;Straub et al., 2019;Smale, 2020). A few studies have shown that grazing and wave exposure may alter and even increase temperature-induced loss of kelp beds (Ling et al., 2009;Rogers-Bennett and Catton, 2019;Butler et al., 2020). Bell et al. (2015) showed that over large biogeographic gradients, decadal climatic patterns (e.g., the North Pacific Gyre Oscillation Index), nitrate availability and wave disturbance largely accounted for giant kelp biomass and had profound spatial differences in effect sizes. Because ecological stressors rarely operate in isolation (Crain et al., 2008;Harvey et al., 2013;McPherson et al., 2021), it is of fundamental importance to identify and test for MHW impacts in concert with other stressors, like low nutrient levels and reduced water clarity. For example, range contraction of giant kelp (Macrocystis pyrifera) near its southern limit in Mexico correlate with cyclical El Niño events because a combination of warm waters and low nutrient levels reduce growth, reproduction and survival (Hernandez-Carmona et al., 2001;Edwards and Hernandez-Carmona, 2005). However, there are few, if any, kelp studies that have formally tested for interactive effects between temperature anomalies, nutrient availability, and water clarity across a latitudinal gradient.
New Zealand, like several other bioregions (Gupta et al., 2020), recently experienced the warmest summer on record associated with the extreme Tasman Sea MHW over the austral summer of 2017-18 when sea surface temperature increased to more than 23 • C (c. 4-5 • C greater than average) around some parts of South Island (Salinger et al., 2019;Salinger et al., 2020). This heatwave caused migrations of northern warmwater fishes, glacial melting, early harvest of summer fruit, and loss of large habitat-forming intertidal southern bull kelp (Durvillaea spp.) (Salinger et al., 2019;Salinger et al., 2020). However, it is unknown if other marine foundation species were affected and if co-occurring stressors modified impacts. In contrast to many other coastlines inhabited by kelp, the southern New Zealand coastline is characterized by high sediment runoff and therefore low water clarity, partly associated with natural conditions (steep topography, soft bedrocks, high precipitation) and partly by extensive land-usage changes, such as converting natural forest to pasture (Ewers et al., 2006;Schiel and Howard-Williams, 2016). It is therefore possible that low water clarity modified ecological impacts from the 2017/18 MHW.
Macrocystis is an iconic and conspicuous, habitat-forming, foundation species inhabiting shallow subtidal rocky reefs along the southern coastline of New Zealand, with a northern range around 41 • S on the southern part of the North Island (Hay, 1990;Nelson, 2020). Macrocystis beds are highly productive and characterized by high biodiversity and extensive surface-floating canopies (Foster and Schiel, 1985;Hepburn et al., 2007;Reed et al., 2016;Mabin et al., 2019;Tait, 2019). Previous research from New Zealand has shown that Macrocystis responds negatively to low light levels and that low water clarity and sedimentation reduce its maximum depth to only 10-15 m (Desmond et al., 2015;Tait, 2019), a compressed range compared to 30-50 m depth distribution in clearer waters in North America (Foster and Schiel, 1985;Graham et al., 2007). Understanding whether recent MHWs and low water clarity in concert affected Macrocystis is critical in understanding and managing effects of sediment loads from land-based sources, as the ocean warms and extreme heat events become more common (Frölicher et al., 2018;Hobday et al., 2018;Oliver et al., 2019).
Here, we quantified the abundances of Macrocystis canopy cover over 5 years, including the most extreme MHW on record, and within 3 zones in each of four regions, covering its latitudinal range on mainland New Zealand. More specifically, we analyzed satellite images over 60 months to determine if sea surface temperature and temperature anomalies affected the cover of surface canopies and if diffuse downwelling attenuation (K d ) and chlorophyll-a (a reasonable proxy for nutrient limitation; Goes et al., 1999Goes et al., , 2000 modified or enhanced temperature effects.

MATERIALS AND METHODS
One consequence of reaching the sea surface is that spectral signatures of Macrocystis are relatively easily detected by satellites Mora-Soto et al., 2020), whereas most other marine macroalgae (except intertidal species) are continuously submerged and their spectral signatures dissipated through the water column. Cavanaugh et al. (2011) presented a thorough example of the requirements for gauging surface canopies over time through remote sensing by satellites; they were able to ground-truth kelp biomass by comparing data from subtidal surveys over different seasons and years to canopy cover using remote sensing products. In New Zealand, unfortunately, there are few estimates of Macrocystis biomass but, nevertheless, remotely quantifying the areal coverage of Macrocystis provides scalable metrics of ecosystem health (Hamilton et al., 2020). One important benefit of passive remote monitoring by satellites is that it removes or reduces the requirement for in situ subtidal surveys, while providing a platform to examine kelp bed size efficiently and objectively over broad geographical ranges and significant oceanographic gradients. However, the limitations and assumptions of using satellite imagery for these purposes must also be recognized.
Passive remote sensing is now widely used Huovinen et al., 2020;Mora-Soto et al., 2020;Hamilton et al., 2020) and the approach generally computes vegetation indices based on measurements of near infrared electro-magnetic radiation. These indices detect vegetation not occluded by overlying water, giving a direct measurement of only the floating portion of macroalgal canopies. By limiting detection to surface canopies, these datasets are unable to integrate the full population dynamics of Macrocystis and give no insight into the presence of subsurface forests. Although these satellite methods cannot observe subsurface forests, surface cover of Macrocystis is a key indicator of the health and vitality of the depth integrated giant kelp populations so that, remotely sensed time series have great potential for identifying the integrated responses of kelp beds to multiple stressors .

Quantifying Macrocystis Cover
Surface canopies of Macrocystis were assessed using Sentinel-2 satellite imagery (resolution = 100 m 2 ) between December 2015 and December 2020 (Copernicus Sentinel-2A data 2015-2020). Four focal regions spanning 6 • of latitude across New Zealand's southern coastline were chosen, relating to key populations (Hay, 1990), and representing the full range of Macrocystis distribution on mainland New Zealand. Within the four focal regions, three zones were selected to encompass as broader range of environmental conditions (e.g., fetch, wave exposure) as possible. Individual beds of key Macrocystis populations were selected to fit within a depth range of 8-15 m. Although data availability for total nitrogen concentrations varied between regions, values ranged between 0.08 and 0.6 mg/L (Dudley et al., 2017). Within each region three zones were selected where Macrocystis was present in 2015. Within each zone, five polygons were allocated (500 × 1,000 m with the long edge of the polygon approximately parallel with the coast) to estimate total coverage of surface canopies. A total of 60 polygons (4 regions × 3 zones × 5 polygons) were thereby used to generate a time series of canopy cover. Kelp cover (m 2 ) was calculated from the number of pixels with detectable vegetation multiplied by the pixel area (100 m 2 ). Data were filtered, masked, and downloaded using Google Earth Engine (Gorelick et al., 2017). Although the timing of satellite capture was not synced to tidal cycles (with a tidal range of ca. 2 m in most of our study region), with potential implications for the visible extent of canopies, remote sensing studies of Nereocystis luetkeana  and Macrocystis pyrifera (Butler et al., 2020) show that variation in coverage was not particularly sensitive to tidal height.
The Normalized Difference Vegetation Index (NDVI) and the Normalized Difference Water Index (NDWI) were calculated from the visible light and NIR bands of the Sentinel-2 constellation. Specifically, the bands B3 (559.8 nm) and B8 (832.8 nm) were used to calculate NDWI and B4 (664.6 nm), and B6 (740.5 nm) was used to calculate NDVI. To differentiate kelp from other features, a series of masking and filtering procedures was performed. First, monthly near cloud-free images were selected by using quality control bands (Sentinal-2 QA band 60). The acceptable percentage of cloud cover was set to 10%. The coastline and offshore islands were masked by an elevation layer to remove all land-based pixels. The thresholds of the NDVI were set >0.01 slightly more conservative than the threshold set for detection of giant kelp by Mora-Soto et al. (2020), but less than thresholds for seagrass (0.2; Calleja et al., 2017), and the threshold of NDWI was set to <0.2 based on detection results of well-studied beds (Tait, 2019). Finally, the monthly cover (when appropriate imagery was available) of Macrocystis was estimated within each of 60 polygons across the four regions from December 2015 to December 2020. See https://leightait.users.earthengine.app/view/ kelpnz for example images, and https://code.earthengine.google. com/62886937d67d2ce5086d76b436c80ce3 for scripts.
Kelp detection results were tested against in situ subtidal densities of Macrocystis (Tait, 2019). In situ subtidal population surveys (Tait, 2019) partially spanned the satellite time-series and allowed near-direct comparisons. We compare the densities of mature Macrocystis measured at eight subtidal sites at 8-11 m depth to satellite estimated coverage in a 100 × 100 m polygon centered at each subtidal site. Comparisons between subtidal densities of Macrocystis and remote estimates revealed a strong positive linear relationship, but also some variation between subtidal densities and the coverage of surface canopies (Supplementary Figure 1).
Although the method does not provide a specific canopy area per pixel, instead assuming 100% canopy coverage within pixels, similar NDVI based vegetation detection methods have been shown to provide an effective proxy for kelp extent and abundance (Cavanaugh et al., 2010;Nijland et al., 2019). This provides a standardized method for identifying the presence and relative extent of kelp beds to identify spatio-temporal trends in relation to key environmental parameters (Butler et al., 2020).
Quantifying Sea Surface Temperature, Light Attenuation, and Chlorophyll-a Sea Surface Temperature Sea surface temperatures (SST) were estimated using the NOAA "Optimum Interpolation Sea Surface Temperature (OISST) V2.1" product (Reynolds and Banzon, 2008). This provides a 1/4 degree global, daily SST estimate from late 1981 to present. Daily SST and temperature anomaly (i.e., the daily OISST minus a 30-year climatological mean) values from December 2015 to December 2020 were extracted from the NOAA OISST product. Daily SST and anomalies were then used to calculate mean monthly SST, mean monthly anomalies and maximum monthly SST for each of the 12 zones. SST means and anomalies were taken adjacent to the areas where kelp abundances were calculated. Data were downloaded using Google Earth Engine (Gorelick et al., 2017).

Chlorophyll-a
Chlorophyll-a concentration (chl-a) was estimated using satellite measurements of ocean color from the Moderate Resolution Imaging Spectrometer on the Aqua satellite (MODIS-Aqua)owned and operated by the US National Aeronautics and Space Administration (NASA, 2018). We used the Quasi-Analytical Algorithm (QAA) algorithm (Lee et al., 2002(Lee et al., , 2009 (488)]. Phytoplankton absorption was converted to an estimate of chl-a using the chl-specific absorption coefficient, aph * (488). The value of aph * (488) can vary seasonally and spatially, relating to different phytoplankton species (with varying cell physiology and pigments), different phytoplankton cell sizes, and the light environment (Kirk, 2011). Here, we used an average of values found for oceanic phytoplankton (Bricaud et al., 1995;Bissett et al., 1997), and measurements in the lower reaches of New Zealand rivers and estuaries. We blended the QAA-chl-a and the MODIS-default chl-a product (NASA, 2018) using a logistic-scaling of bbp(555) (Pinkerton et al., 2018).
Chl-a were estimated from 1 km 2 polygons adjacent to each of our 12 zones and each zone was averaged monthly. Here we used chl-a as a proxy for nutrient availability, particularly nitrogen. While the use of chl-a as a proxy for nutrient concentrations has limitations, it is broadly indicative of nutrient limitation (Goes et al., 1999(Goes et al., , 2000. To do this we compared total nitrogen (including NH 4 + , NO 3 − , NO 2 − ) concentrations to chla concentrations as measured by in situ sampling across much of our study zone (Dudley et al., 2017). These measurements were taken from 30 sites at near-shore locations (between 50 and 500 m from shore) between 2015 and 2017. Depending on location, samples were taken monthly or quarterly. We also investigated the role of temperature on total nitrogen availability (Snyder et al., 2020). The relationship between total nitrogen and chl-a concentrations showed a strong positive correlation (Supplementary Figure 2). However, the negative relationship between temperature and total nitrogen, while significant, had a poor fit (Supplementary Figure 2).

Satellite Estimate of the Downward Attenuation Coefficient of Light (K d )
The diffuse downwelling attenuation coefficient in the Photosynthetically Available Radiation (PAR range, 400-700 nm), K d (m −1 ) was used as our measure of water clarity. Values of K d were estimated from MODIS-Aqua measurements of ocean color, processed to inherent optical properties using the QAA algorithm (Lee et al., 2002(Lee et al., , 2009) following the methodology of Pinkerton et al. (2018). From these IOPs, we estimated the diffuse attenuation coefficient in the PAR range as Lee et al. (2005) and Shi and Wang (2007). The satellite-derived attenuation coefficient was mapped at a nominal resolution of 500 × 500 m and projected to a transverse Mercator grid. The temporal resolution of the product for the study region is 1-2 measurements daily. Values of K d were extracted from the dataset around kelp forests and averaged monthly to provide a dataset with low quantities of missing data (Pinkerton et al., 2018).

Data Analysis
Monthly estimates of kelp coverage from 5 polygons within each zone were averaged and aligned with monthly means of sea surface temperature, water clarity and chl-a estimated at the zone level (12 zones within four regions). Given variation in cloud cover during satellite passes, variation in imagery availability resulted in variable sample numbers across regions (North = 160, North-East = 135, South-East = 169, South = 132).
The effects of monthly maximum SST, temperature anomaly, water clarity (as defined by the light attenuation coefficient K d ) and chl-a concentration (as a surrogate for nutrient availability) on Macrocystis cover were analyzed with Generalized Additive Models (GAMs) using the "R" package "mgcv." Furthermore, GAMs were used separately to assess the potential for linear and non-linear covariance between physical parameters (temperature anomalies, light attenuation, and chl-a concentration). Assumptions of normality (Q-Q plot), homogeneity of variance (Levene's Test), as well as "concurvity" for general additive model analysis (an estimate of redundancy among explanatory variables) were checked for models. GAM models were fitted with "cr" (cubic regression) splines, using a k-value of 6 (i.e., the number of "knots" denoting the complexity of the non-linear fit), and the distribution family "tweedie." Selection procedures were implemented to penalize and remove factors with poor explanatory power. The final model included mean monthly temperature, K d , temperature anomaly, chl-a, and the two-way interaction between water clarity and temperature anomalies. Furthermore, the categorical variable "zone" was included in the GAM model.

RESULTS
Mean SST anomalies between December 2015 and December 2020 showed strong warming trends along much of the coastline of southern New Zealand (Figure 1A). The normally cooler southern regions had the greatest mean anomalies (+1.5-2 • C warmer than average) between 2015 and 2020. Conversely the normally warmer northern regions experienced some of the coolest mean temperature anomalies (−0.5-0 • C cooler than average) across the whole region ( Figure 1A). The influence of mean monthly temperature anomalies on the mean monthly cover of Macrocystis varied across regions, with minor warm anomalies (c. +1 • C) reducing the kelp cover toward the northern limits of its range, whereas at the southern-most region (4. South), large anomalies (c. +3 • C) had a negative but less dramatic effect on coverage ( Figure 1B). In the central regions (2. North-East, 3. South-East) warm anomalies associated with summer MHWs caused notable reductions in surface canopies ( Figure 1B). The central region (2. North-East) had the most extreme warm anomalies (4.7 • C above average for December 2017, during the severe MHW), and the biggest range of anomalies (Supplementary Figure 3). The southern two regions had no monthly mean temperatures below average for the 5-year period (Supplementary Figure 2).
At the northern extent of Macrocystis, SST averaged 16 • C over the 5-year period (December 2015 till December 2020), but at the southern boundaries of giant kelp the average dropped to c. 13.5 • C (Figure 2A). The mean and full range of temperatures experienced (as shown by the span of the x-axis, Figure 2B) within regions varied over the 5 years, the North region averaged c. 16 • C and spanned 9 • C, North-East averaged c. 15 • C and spanned 11.5 • C, South-East averaged had c. 14.0 • C and spanned a 9.5 • C range, and the South region averaged c. 13.5 • C and had the narrowest thermal range (7.5 • C) ( Figure 2B). Declining cover of Macrocystis was observed in all regions, but by far the biggest response was in the North-East and South-East regions where temperatures above c. 18 • C caused large reductions of Macrocystis canopy cover ( Figure 2B).
The effect of water clarity on kelp coverage varied across regions, as did the range of light attenuation observed. Macrocystis populations in the North (1) experienced uniformly higher water clarity (i.e., K d ), with no detectable response in kelp cover relating to it as a single factor ( Figure 3A). However, the North-East and South-East regions (2-3) experienced the greatest range of water clarity, which related to a loss of canopy cover. This relationship was clear only for the North-East region (2), however, as the apparent negative relationship in the South-East region was dependent on only one reading of high K d and did not depict a reliable trend (Figure 3B). Conversely kelp coverage of the South region (4) was, somewhat anomalously, positively related to increasing K d (Figure 3B). This region had a positive relationship between kelp coverage and potential nutrient availability (i.e., chl-a; Figure 3B), but this was not statistically significant.
The time series of kelp coverage averaged across all regions and zones showed a dip in kelp coverage that was coincident with the severe MHW in 2017-2018 (Figures 4A,B and  Supplementary Figure 3). Water clarity and chl-a showed no discernible temporal trend (Figures 4C,D), although there was a slight depression of chl-a during the 2017-2018 MHW. To assess covariance among predictor variables (chl-a, K d , temperature anomalies) we assessed the strength of linear and non-linear correlation using a separate GAM model excluding the influence of predictors on kelp coverage. Covariance analyses of these variables showed strongest correlations between chl-a and K d , and temperature anomalies and chl-a ( Table 1). There was no significant correlation between temperature anomalies and water clarity ( Table 1). Across all regions, the coverage of kelp showed strong trends with temperature, and K d (Figure 5). Temperature anomalies had a peak positive effect on kelp coverage between +1 and +2 • C, and a pronounced negative effect beyond +2 • C ( Figure 5A). Maximum monthly temperatures between 9 and 18 • C were associated with stable coverage of kelp, but coverage declined steeply beyond 18 • C (Figure 5B). Kelp coverage was highly variable at the low ends of K d (clearer water) and declined beyond K d values of 0.3 ( Figure 5C). Increasing potential nutrient availability (as indicated by chl-a) had no discernible overall effect on kelp cover ( Figure 5D). While interactions between chl-a and other parameters were tested, none contributed to deviance explained and were dropped from the final model. General additive models (GAM) indicated that temperature anomalies, maximum monthly temperature, and light attenuation contributed most to the cover of Macrocystis (Table 2 and Supplementary Figure 2). Concurvity, or collinearity of variables showed that model parameters were between 0.5 and 0.17 (values close to one indicate redundancy). The interactive model also revealed a significant two-way interaction between temperature anomalies and K d ( Table 2 and Supplementary Figure 4). The interactive effect between temperature anomalies and K d was shown in a two-way interactive model, where rising temperature anomalies and falling water clarity combined to strengthen the decline in kelp coverage (Figure 6). In clear water (i.e., low K d ) temperature anomalies had little effect on kelp cover between −2 and +1 • C, but beyond +1 • C kelp coverage declined more rapidly. This effect was accentuated in low water clarity (i.e., high K d ).

DISCUSSION
Giant kelp, Macrocystis pyrifera, is one of earth's fastest growing photosynthetic organisms, and is a key contributor to carbon fixation and habitat provision for temperate marine ecosystems across a large extent of the world's temperate coastlines (Schiel and Foster, 2015;Miller et al., 2018). Recent studies have shown that Macrocystis and other large kelps have had retractions or reductions in their distribution in the northern hemisphere (Arafeh-Dalmau et al., 2019;Rogers-Bennett and Catton, 2019) in response to environmental fluctuations and combined physical and trophic interactions (Ling et al., 2009(Ling et al., , 2015Bell et al., 2015). In southern Australia, Macrocystis has experienced massive retractions, and is nearing functional extinction in some regions in response to oceanographic shifts affecting nitrogen concentrations, temperature, and larval delivery of a key sea urchin grazer (Johnson et al., 2011;Mabin et al., 2019;Butler et al., 2020). The results of our study suggest that if current seawater warming trends continue, they could affect the cover of Macrocystis across much of its current range in southern New Zealand. Although temperature was the major driver we identified, it is also evident that declining water clarity has negative effects on giant kelp. We also showed that the negative effect from high water temperature is exacerbated by low water clarity. Unlike other studies (Dean and Jacobsen, 1986;Mabin et al., 2019), however, the correlation between MHW events and depressed chlorophyll-a (as a proxy for nutrient limitation) was not evident. The southernmost region studied (−46 • latitude) showed some resilience to warming events but those just to the north (−45 • latitude) were more affected. This may have been due to the intensity of the MHWs in those regions, which highlights that MHWs do not have a uniform and constant effect along a heterogeneous coastline. The steep decline with increasing temperatures in the central regions may be due to potentially greater susceptibility of cold-adapted populations. At the northern extent of Macrocystis distribution (-41 • S) temperatures during the 5-year study period were colder on average (compared to the 38-year time series of sea surface temperature; Reynolds and Banzon, 2008), but even small positive anomalies were associated with reductions in cover. This is surprising given that the extreme warm temperatures experienced in New Zealand are 3-5 • C cooler than those experienced by Macrocystis forests in California (2014California ( -2016Cavanaugh et al., 2019). Some evidence suggests that higher latitude populations are more vulnerable to reproductive failures under temperature stress (Hollarsmith et al., 2020). Additional stressors may also be responsible for this massive disparity in temperature thresholds. Of the parameters tested (expect temperature) at a larger scale over all regions, only declining water clarity (increasing K d ) was associated with loss of kelp canopies. In situ studies within the South-East region (3) showed declines in all life history stages of Macrocystis at sites with locally compromised light availability (Tait, 2019). Although the surface canopies detected during this study are less immediately threatened by low water clarity, the ability of juvenile gametophytes and sporophytes to grow and reach the surface will be compromised at light limited locations, or during prolonged low light events. Our study shows that the consequences of warmer temperatures are greater in conditions of compromised light where greater metabolic demand (Mabin et al., 2019) is not met by sufficient light delivery to the benthos (Tait, 2019).
Unlike other studies (e.g., Snyder et al., 2020), we show poor correlation between SST and nitrogen concentration. There was, however, good agreement between in situ measurements of total nitrogen and chl-a concentration for sites covering much of the study region. We attribute this lack of alignment between this study and other research (Snyder et al., 2020) to the high Overall model had an adjusted r 2 = 0.45, deviance explained = 60%, n = 596. Bold and italics refer to statistically significant results.
riverine influence along the southern New Zealand coast and large-scale land-use changes (Schiel and Howard-Williams, 2016;Stevens et al., 2021) and the coastal focus of this study. While the chl-a concentration (Pinkerton et al., 2018) and water clarity algorithms used in this study have been developed to better account for the often turbid waters of southern New Zealand, these case-2 coastal waters present challenges for determining metrics of both light availability (K d ) and phytoplankton biomass. However, the utility of remote sensing to retrospectively assess the effects of multiple physical and biotic parameters on kelp beds has great utility to better understand the interaction between local and global scale stressors. Studies in Central and Southern California have shown highly correlated occurrences of low nitrogen and high temperature (Snyder et al., 2020), while we show some correlation between chl-a (as a reasonable proxy for nutrient availability) and temperature, but not K d and temperature, suggesting that MHW events in New Zealand possibly influence nutrient limitation but not water clarity. There was, however, little indication that nitrogen levels frequently fall below thresholds for nutrient limitation of Macrocystis (Dean and Jacobsen, 1986). Laboratory based assessments of multiple photo-physiological metrics across individual stressor variables (temperature, light availability, and nitrogen availability) clearly show negative impacts of warm temperatures, low light availability, and low nutrient concentrations (Mabin et al., 2019). However, the growth of juvenile sporophytes was not affected greatly by nitrogen concentrations (Mabin et al., 2019). Specific combinations of conditions experienced by real-world populations can significantly affect the outcomes to macroalgae. For example, nutrient history prior to heat-wave events has significant bearing on the outcomes for Macrocystis populations (Fernández et al., 2021), suggesting that the timing and length of heat-wave events or conditions of limited nutrients have significant effects on kelp bed persistence. Our results indicate that chl-a is a reasonable proxy for total nitrogen concentrations, but we present little evidence of major declines in Macrocystis associated with periods of low chl-a (as detected remotely). Despite the fundamental relationship between warm temperatures and low nutrient concentrations (Snyder et al., 2020), the combined evidence of this study and others (Cavanaugh et al., 2019;Butler et al., 2020) is that these kelp forests are highly susceptible to MHWs. We show that reductions in water clarity present a major FIGURE 6 | Interactive effects of temperature anomalies and light attenuation on surface canopy cover of Macrocystis pyrifera across New Zealand. Plot shows the mean response. Contours show the mean response of kelp coverage to each variable, extrapolated 25% beyond available data. Contour plots fitted by general additive models (GAM).
challenge to the resilience of Macrocystis at elevated temperatures. Few studies have used remotely sensed metrics of water clarity to assess the consequences to remotely sensed kelp beds over large spatial scales, and our research presents compelling evidence that regions of low water clarity increase the susceptibility of kelp beds to heat waves.
The use of satellite remote monitoring for tracking populations of Macrocystis is limited to detecting the surface canopies as populations at the benthos cannot be sighted but nevertheless enables kelp forest health and abundance to be monitored remotely over large areas. There is, however, a lag between environmental change and Macrocystis response, so that legacy effects of significant oceanographic and climatic events require further evidence to quantitatively tie-down the independent effects of multiple variables. For example, numerous studies point to the key role of significant wave height in affecting kelp bed dynamics Reed et al., 2011;Bell et al., 2015;Castorani et al., 2018). Synthesis of long-term observations using satellite remote sensing revealed that recovery of Macrocystis canopies following winter storms was highly correlated to spring sea surface temperatures (SST) . Remote tracking of Macrocystis populations has enabled analyses of broad scale processes relating to oceanographic trends Hamilton et al., 2020), and these methods have equal promise for better understanding the point-source and diffuse influences of land-use changes.
New Zealand has experienced some of the most intense marine heat waves on record in the past 5 years (Salinger et al., 2019(Salinger et al., , 2020 causing localized losses of southern bull kelp . Rapid rates of land-use change associated with agriculture and urbanization have greatly altered the land-water interface globally, including in New Zealand, where rates of sedimentation have also increased (Goff, 1997). The proximity of Macrocystis forests in the nearshore zone to sources of sediments greatly affects the demography of populations of large brown algae that are the facilitators of diversity and energy flow in nearshore waters. Sediments may prevent attachment of algal propagules to benthic surfaces Schiel, 2003: Schiel et al., 2006), smother those that have settled (Schiel and Gunn, 2019), and result in a poor conversion from juvenile Macrocystis sporophytes to adult plants reaching the surface (Tait, 2019). The consequences of poor landmanagement and high levels of sediment runoff must urgently be addressed in the face of increasing occurrences of MHWs (Schiel and Howard-Williams, 2016).
Declines in Macrocystis have been ascribed to anomalous warm events, but other mechanisms for kelp loss include altered trophic interactions and destructive grazing by urchins (Ling et al., 2009;Rogers-Bennett and Catton, 2019;Butler et al., 2020) and nutrient limitation (Dean and Jacobsen, 1986;Hernandez-Carmona et al., 2001;Edwards and Hernandez-Carmona, 2005;McPherson et al., 2021). Many regions, however, such as southern New Zealand, have relatively few instances of overgrazing by urchins (Schiel and Foster, 2015). At this stage, it is unknown how many other regions are experiencing consequences of anthropogenically altered sediment regimes on the abundance and distribution of kelp forests. This is likely to be exacerbated with global warming, as increasing loss of glaciers (Smith et al., 1999) threatens to increase transport of terrigenous sediments to the oceans. In temperate regions where coastal marine ecosystems meet managed terrestrial environments, appropriate interventions to reduce the rate of terrigenous sediment inputs may greatly improve outcomes for giant kelp forests and enhance their resilience to threats associated with MHWs and other global or regional stressors. Failure to act in supporting temperate kelp forest ecosystems has major implications to the operation of global carbon sink pathways (Filbee-Dexterand Wernberg, 2020), and overall ecosystem stability.

AUTHOR CONTRIBUTIONS
LT conceived the research, prepared and analyzed data, and drafted the manuscript. FT and MP prepared and analyzed data and to the draft manuscript. MT and DS contributed to the research concept and to writing the manuscript. All authors contributed to the article and approved the submitted version.