Multi-decadal decline in cover of giant kelp Macrocystis pyrifera at the southern limit of its Australian range

: Knowledge of long-term and multi-scale trends in ecological systems is a vital component in understanding their dynamics. We used Landsat satellite imagery to develop the first long-term (1986− 2015) data set describing the cover of dense surface canopies of giant kelp Macrocystis pyrifera around the entire coastline of Tasmania, Australia, and assessed the extent to which potential environmental drivers explain the dynamics of surface canopies at multiple spatial and temporal scales. Broad-scale temporal patterns in canopy cover are correlated with El Niño−Southern Oscillation events, while regional patterns are related to sea surface temperature and nutrient regimes are associated with the East Australian Current. Regression models developed to predict the presence or absence of giant kelp canopy emphasise the importance of sea surface temperature in these systems. Long-term decline in canopy cover is clearly evident in most regions, and in light of increasing thermal stress associated with a changing ocean climate, this raises concern for the future of this species as a major habitat-forming kelp in Australia and some other regions worldwide. Given that Tasmania represents the stronghold of the range of this species in Australia, but is a geographic trap in that there is no suitable habitat for M. pyrifera to the south, our findings support the Federal listing of giant kelp communities in Australia as an endangered marine community type.


INTRODUCTION
Many marine ecosystems are under threat as a consequence of global warming (Hughes et al. 2003, Wernberg et al. 2016. Oceans are becoming warmer and more acidic, shifts in ocean circulation are altering temperature and nutrient regimes, and in many locations, anthropogenically associated disturbance events are increasing in frequency and severity, so that marine ecosystems are increasingly subject to multiple environmental stressors. These physical changes to the environment affect biotic communities, including in coastal regions where eco systems are particularly subject to growing anthropogenic pressure from urbanisation, development and exploitation (Vásquez 2008, Strain et al. 2014). Under-ABSTRACT: Knowledge of long-term and multi-scale trends in ecological systems is a vital component in understanding their dynamics. We used Landsat satellite imagery to develop the first long-term (1986− 2015) data set describing the cover of dense surface canopies of giant kelp Macrocystis pyrifera around the entire coastline of Tasmania, Australia, and assessed the extent to which potential environmental drivers explain the dynamics of surface canopies at multiple spatial and temporal scales. Broad-scale temporal patterns in canopy cover are correlated with El Niño−Southern Oscillation events, while regional patterns are related to sea surface temperature and nutrient regimes are associated with the East Australian Current. Regression models developed to predict the presence or absence of giant kelp canopy emphasise the importance of sea surface temperature in these systems. Long-term decline in canopy cover is clearly evident in most regions, and in light of increasing thermal stress associated with a changing ocean climate, this raises concern for the future of this species as a major habitat-forming kelp in Australia and some other regions worldwide. Given that Tasmania represents the stronghold of the range of this species in Australia, but is a geographic trap in that there is no suitable habitat for M. pyrifera to the south, our findings support the Federal listing of giant kelp communities in Australia as an endangered marine community type. standing how these changes manifest, however, is difficult. Species' responses are the re sult of complex biotic and abiotic interactions, and long-term datasets are required to separate anthropogenically driven changes from natural variability (Reed et al. 2016. Data collection in the marine environment is costly and often labour intensive, so despite recognition of the value of long-term data sets, few studies are able to assess long-term change.
Among coastal ecosystems, giant kelp (Macrocystis pyrifera) forests are a key and iconic habitat that dominate many nearshore rocky coastlines in temperate and cold-water regions worldwide (Graham et al. 2007). Growing from the seafloor and extending as much as 40 m to the surface, these dense underwater forests form closed to semi-closed canopies that alter the light environment, reduce water flow, alter sedimentation rates and provide an important source of fixed carbon both within the kelp forest and to surrounding communities (Steneck et al. 2002, Layton et al. 2019. Their 3-dimensional structure provides habitat for a diverse range of fish and invertebrate species, including as recruitment grounds for commercially and non-commercially exploited species, and constitutes one of the most biodiverse marine systems in the world (Schiel & Foster 2015). Many of these species are directly associated with the forest itself, and loss of the canopy can reduce the trophic structure of the ecosystem (Graham 2004, O'Connor & Anderson 2010 and have significant impacts on biodiversity. The success of giant kelp populations is in part due to their ecological plasticity. Given an adequate spore supply, suitable rocky substrata and a range of environmental conditions, giant kelp can rapidly establish, grow and become reproductive in less than a year ). These qualities have enabled it to colonise and adapt to coastal habitats across the globe, an ability that is further enhanced by its capacity to float and disperse via rafting (Rothäusler et al. 2009) and its capacity for genetic adaptation (Camus et al. 2018).
Reproduction, recruitment and growth in giant kelp are controlled by physical ocean properties, including nutrient concentrations and temperature (Graham et al. 2007). The limited capacity of M. pyrifera to store nitrate, a nutrient essential for photosynthesis, means that populations will respond rapidly to changes in nitrate and ammonia concentrations (Kop czak 1994). Microscopic and juvenile sporophytes demonstrate a similar response, and recruitment occurs only when favourable conditions of nitrogen availability and light coincide (Dean & Jacobsen 1984, Carney & Edwards 2010. Temperature is almost always negatively correlated with nitrate (although see Fram et al. 2008, Brzezinksi et al. 2013, and giant kelp populations do not typically survive or establish in waters > 20°C (Cavanaugh et al. 2019). A growing body of literature is also reporting that early life stages of giant kelp are vulnerable to the effects of temperature even when nutrients are abundant (e.g. Hollarsmith et al. 2020).
In areas where conditions are otherwise suitable for growth, wave-driven disturbance is a primary control of mortality in M. pyrifera sporophytes . Because of their large size, giant kelps experience high drag forces which make canopies and even whole plants vulnerable to dislodgement during periods of high water turbulence (Seymour et al. 1989). Prolonged periods of low nutrients and/or high temperatures, such as those associated with El Niño−Southern Oscillation (ENSO) events, can also lead to frond deterioration, senescence or sporophyte mortality (Dayton et al. 1999, Edwards 2004. Population dynamics in giant kelp thus reflect the interplay between resource availability (nutrients, temperature) and physical disturbance (wave activity), as determined by seasonal and inter-annual fluctuations in ocean circulation, as well as latitudinal gradients and variability driven by local seafloor topography.
A number of studies have documented local and regional changes in the distribution and extent of giant kelp populations and their relationship with inter-annual climate fluctuations. These studies consistently record severe decline and elimination of giant kelp forests during El Niño periods in the eastern Pacific when coastal upwelling is depressed, but recovery is typically observed within <1 to 4 yr. For example, Edwards (2004) examined the effect of the 1997− 1998 El Niño on populations of giant kelp over its entire geographic range in the northeast Pacific, and found that while the stormy, warm and nutrientpoor conditions associated with these events resulted in the near-complete loss of giant kelp in the southern half of its range, recovery occurred within 6 mo to 2 yr. On the Mejillones Peninsula in Chile, disturbances associated with this same 1997−1998 El Niño led to the loss and subsequent recovery of M. pyrifera (morph integri folia) populations within 4 yr (Vega et al. 2005), while in Mexico (Baja California) and Chile, M. pyrifera forest canopy cover recovered to former levels within 3 yr of significant losses (Edwards & Hernández-Carmona 2005. Despite the importance of assessing short-term dynamics of M. pyrifera in relation to climate fluctuations, few studies have been able to assess multidecadal trends, in part due to a lack of long-term data. The release of free satellite imagery by the US Geological Survey (USGS) in 2008 has made available a data source of unprecedented spatial and temporal resolution, dating back to 1982. A prominent use of this data archive was the work of Cavanaugh et al. (2011), who used a 25 yr time-series of images to investigate biomass dynamics in giant kelp populations in the Santa Barbara Channel (California, USA). They revealed significant correlations between biomass and annual and inter-annual climate signals at regional scales, but little overall change in biomass over the entire time period. These results are echoed in ongoing research using the same Landsat archive in this and the broader Californian region (e.g. Bell et al. 2020), but whether these longterm trends and links between biomass and climate signals are similar in other locations around the world is unknown.
The coastal waters of the state of Tasmania represent the most extensive giant kelp habitat in Australia. Large and dense forests have historically covered significant areas of the nearshore reef habitat and formed a key and characteristic component of coastal ecosystems, particularly in eastern Tasmania, but significant declines have been recorded in eastern Tasmania in recent decades. In one of the first long-term assessments of abundance in the state, Johnson et al. (2011) reviewed aerial photographs of 7 sites along the east coast, and reported up to 98% declines in canopy cover between 1946 and 2007. This was the first study to quantify the loss of giant kelp in Tasmania, and prompted the 2012 listing of the giant kelp community through the Federal Environment Protection and Biodiversity Conservation Act as an 'endangered marine community type'. Given this listing, the first ever of a marine community in Australia, a state-wide assessment of trends in giant kelp canopy cover is of particular importance to address conservation concerns.
In the present study, we used satellite imagery to quantify the extent of giant kelp surface canopies for the entire coastline of Tasmania over 30 yr, from 1986− 2015. We relate the spatial and temporal variability in surface cover of M. pyrifera to ocean environmental conditions known to affect giant kelp populations (temperature, exposure), and quantify whether these factors influence the presence of M. pyrifera canopies in this region. Our primary aim was to assess space−time variability in this ecologically important and iconic kelp habitat, filling a gap in knowledge for the region and contributing to a growing global dataset on the ecology of giant kelp.

Study area
Tasmania is an island state located off the southeast corner of Australia (Fig. 1). The waters surrounding Tasmania offer a highly variable environment in terms of oceanic regimes, with different areas of the coastline subject to different ocean current and exposure conditions that vary according to seasonal and inter-annual ENSO cycles. The most dominant oceanic currents to impact Tasmania are Australia's east and west continental boundary currents, the East Australian Current (EAC) and the Zeehan Current (ZC; the southern-most extension of the Leeuwin Current). Both currents are characterised by warmer, high-salinity but nutrient-poor water that flows southwards along the east (EAC) and west (ZC) coasts (Ridgway & Condie 2004, Oliver & Holbrook 2018. In Tasmania, the EAC is an erratic eddy-driven flow that is strongest in summer, while the ZC demonstrates a more consistent flow that peaks in winter (Oliver & Holbrook 2018). In transitional stages where neither current dominates, the influence of incursions of the Antarctic Circumpolar Current (ACC) are recognised by its cold and relatively nutrient-rich signature (Bucha nan et al. 2014).
Seasonal cycles in temperature and nutrients are superimposed onto inter-annual cycles corresponding to the ENSO. ENSO describes changes to the circulation of the South Pacific Gyre, which cycles through El Niño and La Niña events every 8−12 yr. La Niña events result in a strengthening of the EAC and ZC and subsequent warming of coastal waters off Tasmania, while El Niño events are characterised by weakened currents and cooler waters (Feng et al. 2003, Oliver & Holbrook 2018.
Exposure conditions also vary across the state. The west and south coasts are subject to frequent and intense weather systems originating in the Southern Ocean that produce strong winds and large waves (Hill et al. 2010). The north coast is offered some protection from within Bass Strait, but experiences strong tidal currents and wind-generated waves (Porter-Smith et al. 2004). The east coast is subject to variable weather conditions (Porter-Smith et al. 2004, Hill et al. 2010, with occasional storm events generated from low pressure systems to the east.

Quantifying the extent of and trend in surface canopies of Macrocystis pyrifera
2.2.1. Satellite estimation of M. pyrifera canopy Giant kelp can form large and dense floating canopies at the sea surface. These canopies are photo synthetic and have spectral properties very similar to those of terrestrial vegetation, namely a low re flect ance in the visible range and a high reflectance in the near infrared (NIR). Water absorbs almost all of the incoming NIR wavelengths, and it is this contrast that enables identification of M. pyrifera canopies in satellite imagery.
The presence of surface canopies of M. pyrifera was estimated from Landsat TM, ETM+ and OLI Sur-face Reflectance Product imagery. Landsat satellites have been in operation almost continuously since 1984, and have acquired multispectral imagery with a 30 m spatial resolution and a 16 d repeat cycle across the globe. This imagery includes bands in both the visible and NIR ranges. Images were sourced from the USGS (www. earthexplorer. usgs. gov), where an atmospheric correction procedure had already been applied (for details see USGS 2016USGS , 2017. Images were classified based on the normalised difference vegetation index (NDVI). The NDVI is commonly used to detect terrestrial vegetation and has more recently been shown to provide an effective proxy for kelp extent and abundance (e.g. Cavanaugh et al. 2010, Nijland et al. 2019. The NDVI describes the ratio of NIR to red (R) reflectance of a given pixel as: (1) A threshold NDVI value was selected such that when it was ex ceeded in a given pixel, that pixel was classified as containing 100% giant kelp canopy floating at the surface. Although this method does not provide a quantification of the specific canopy area within a given pixel, it provides identification of the presence of surface canopies within the area and an adequate means to assess spatio-temporal trends in forest habitat. The threshold was chosen based on how well it delineated canopies that were observed in the false-colour satellite image, across multiple images. Although numerous canopy-forming kelps exist in Tasmania (e.g. Ecklonia radiata, Phyllospora comosa, Durvillaea potatorum), giant kelp is the only surface canopy-forming species in waters at 5−40 m depths, making identification of this species unambiguous.
The slight differences in the band designations between Landsat 4/5 TM, Landsat 7 ETM+ and Landsat 8 OLI has minimal effect on the NDVI (Ke et al. 2015, T. Bell pers. comm.). Strong positive relationships (r 2 > 0.994) in kelp cover fractions derived from the different sensors have also been demonstrated using the MESMA technique . Considering that the focus of this work is on detecting presence or NDVI NIR R NIR R = − + absence of surface canopy cover, and not a quantitative analysis of NDVI, no correction to the NDVI was required.
In the marine environment, the reflectance signal, and thus the NDVI value, can be affected by sun glint, breaking surface waves, phytoplankton blooms, dissolved organic matter and sediment run-off. In such cases where this leads to noise being introduced into the classification product, manual interpretation is necessary, and in our case was validated against the presence of canopy cover as ob served in the falsecolour satellite image. Although a reasonably simple means of validation, the use of NDVI for mapping kelp surface canopies has been applied and validated in other regions worldwide and is considered an effective means of mapping trends in canopy extent through time , Nijland et al. 2019.
The region of each image that was examined for M. pyrifera canopies was between 30 and 1000 m horizontal distance from the shoreline. This defined the area within which it was likely for M. pyrifera to occur. The inshore limit of 30 m was chosen because although M. pyrifera often occurs within 30 m of the coast, there are other algal species that form canopies in the intertidal and shallow subtidal zone (e.g. D. potatorum) that might potentially lead to misinterpretation, particularly if the image was taken at low tide. The offshore limit of 1000 m defines the outer ex tent of probable giant kelp habitat, as be yond this, water depths are typically > 30 m where surface light penetration is inadequate for growth in most cases (Schiel & Foster 2015). The zone was delineated using a mask generated in ArcGIS (ESRI 2014).
The 4 km wide (and shallow) channel between Acteon Island and mainland Tasmania is well known for dense surface canopies of giant kelp; however, as this region is situated outside of the delineated zone (i.e. >1000 km offshore), the mask was extended in this location to ensure that this important area for giant kelp forests was included in the analysis (see Fig. 1).

Data collection.
The extent of and trend in surface canopies of giant kelp were assessed using (1) an inventory of canopy cover across the entire coastline at a coarse (3 yr) temporal resolution, and (2) a fine-scale (2 wk to 12 mo) time-series of canopy cover collated for 24 individual sites distributed around the state. For (1), the 3 yr temporal resolution was necessary to ensure complete cloud-free cover-age of the coastline during the period of maximum canopy development between July and December.
2.2.2.2. State-wide extent of canopy cover. Eight maps at 3 yr intervals, commencing in 1987, were generated to show the extent of dense surface canopies around the entire state through time (herein referred to as the state-wide analysis). Only images taken between 1 July and 30 December were considered, as this represents the period of maximum canopy growth in this region. Where more than 1 cloud-free image existed for an area, preference was given to images with high quality ratings that were as close to September (the spring growth maximum) as possible. The 90 Landsat scenes selected for this analysis and associated metadata are outlined in Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ m653 p001_ supp. pdf.

Trends in abundance of canopy cover.
Trends in the abundance of giant kelp canopy cover were investigated through building a time-series at a much finer temporal resolution. We chose 24 spatially defined sites around the state based on the presence of M. pyrifera canopy in the state-wide ana lysis (Fig. 1). Due to the small total area of kelp canopy observed on the northern Bass Strait coast in the state-wide analysis, this region (corresponding to the Otway, Boags and Flinders bioregions; see Fig. 1) was not investigated at this finer temporal and spatial scale. The final number of sites across the state reflected the number of localised areas where M. pyrifera canopies were observed with reasonable recurrence through time in the state-wide analysis. All cloud-free images between July and December were considered. This fine-scale analysis produced a total of 1068 observations of surface canopy cover across the 24 sites from 1986−2015. Landsat scenes selected for this analysis and associated metadata are outlined in Table S2.
The selected images were classified and analysed using ENVI version 5.2 (Exelis Visual Information Solutions) and ArcGIS (ESRI 2014) softwares. Once classified, the area of M. pyrifera canopy was calculated, and canopy locations were sorted into the 7 bioregions of Tasmania, de fined by the Integrated Marine and Coastal Regionalisation for Australia (IMCRA Technical Group 1998). The IMCRA bioregionalisation is based on characteristics of marine community assemblages that correlate with exposure and ocean regimes.
2.2.2.4. Assessing sensitivity to tide height. Tide height can affect the extent of canopy present at the ocean surface (Britton-Simmons et al. 2008). We examined effects of tide on the area of kelp canopy ob -served by comparing pairs of images taken at the same site that were close in time (16 d apart, the minimum time it takes for Landsat satellites to pass over the same location). Hourly tide height estimates were obtained for each individual site in the fine-scale ana lysis for the period of the study from the FES2012 tide model. Pearson's correlation coefficient was used to validate model-derived tide heights against observations obtained by the Bureau of Meteorology at Spring Bay on the east coast of Tasmania over the same time period (r 271751 = 0.99, p < 0.001). FES 2012 tide height was extracted for each observation by rounding the time of Landsat image capture to the nearest hour to match tide height predictions. Heights were categorised into either low, mid or high tide by first calculating the difference (d) in maximum to minimum tide height at each site. Low tide was classified as the lowest third of tide heights (min + ¹ -3 d), mid tides as the middle third (values between min + ¹ -3 d and max -¹ -3 d), and high tide as the highest third (max -¹ -3 d). Paired images were divided into 2 groups based on whether there was a difference in tide height category be tween the initial observation and the second image 16 d later, or whether the tide height category was the same for the first and second image. The percentage change in canopy area from (1) lower tide to higher tide, or (2) earliest observation to latest observation was calculated. Students t-tests were used to test (1) whether the percentage change in canopy area (from highest tide category to lowest tide category) was significantly different from zero when there was an observed change in tide height category, and (2) whether the average percentage change in canopy area differed significantly between the 2 groups (i.e. when there was a change in tide height category, and when there was not a change in tide height category). Table S2 lists tide heights, estimated as above, at the time of image acquisition.

Sea surface temperature (SST)
SST is important in influencing survival and growth of giant kelp in its own right, but in many cases can also be considered as a proxy for nutrient availability, particularly nitrate, which also affects the growth and cover of surface canopies of giant kelp. Daily time-series of SST for each site were extracted from the NOAA National Climate Data Centre Optimally Interpolated AVHRR product for the period 1981−2016 at a 1/4° grid resolution (Optimum Interpolation Sea Surface Temperature, www. ncdc. noaa.gov/oisst). Stobart et al. (2016) performed a comparison of in situ and satellite SST data in Tasmania, and concluded that satellite-derived measurements provided an adequate estimation of subsurface conditions.
To assess spatial and temporal variability in SST around the state, it was necessary to remove seasonal variation from each time-series to prevent it from confounding any potential difference between bioregions. Each time-series was de-trended by fitting a linear model to each bioregion, and subtracting the slope of the regression from the original data. A sine curve was fitted to the de-trended data, and the resulting function was subtracted from the original data, leaving a daily time-series of mean SST for each bioregion with the long-term trend but without seasonal variation. Long-term trends in SST were asses sed using ANCOVA, with time as a covariate.

ENSO
The Southern Oscillation Index (SOI) gives a measure of the strength and severity of ENSO in Australia. Strong ENSO events alter nutrient, temperature and exposure regimes over periods of 4−7 yr, and each of these variables is important in controlling the biomass dynamics of giant kelp across the globe (Graham et al. 2007). Although there is inherent correlation between SST and ENSO, in southeast Australia, the relationship is typically weak (r ~ 0.3, although statistically significant at the 95% confidence interval; Holbrook & Bindoff 1997).
A monthly time-series of historic SOI was sourced from the Bureau of Meteorology (www.bom.gov.au/ climate/ current/ soihtm1.shtml). The Bureau calculates the SOI as the standardised anomaly of the mean sea level pressure difference between Tahiti and Darwin. The SOI value is calculated on a monthly basis for the whole Australian continent, including Tasmania, and this value was used directly in this analysis.

Exposure
Exposure was estimated at each site using a wave energy model (Keane et al. 2019). Briefly, points were generated at 200 m intervals along any stretch of coastline that was within 200 m of observed giant kelp canopy. For each point, 48 radiating fetchlines were generated, spaced at 7.5°, and extended to a maximum length of 650 km unless intersected by land. Data on swell direction, significant wave height and wave period were then extracted from a 30 yr high-resolution wave hindcast model available from the Bureau of Meteorology (www.bom.gov.au/ climate/ data-services/ocean-data.shtml). This model produced an hourly time-series of wave parameters in a 4 arc-minute gridded output, where wave statistics included both wind and swell components of the wave. Bathymetry was also taken into account in the calculations. For each point and each hourly time step, the model-derived wave parameters were extracted from the bilinear interpolation of the 4 grid cells nearest to that point. The direction of the swell was then rounded to the same set of 48 × 7.5° lines as used for the fetchlines. The normalised fetchline length (normalised to 1) was taken from the fetchline matching the rounded swell direction and multiplied by the derived wave energy (WE = 0.57 × hs 2 × tp, where hs = significant wave height [m] and tp = wave period [s]) to give a wave energy index (WEI). We used this final hourly WEI to generate a mean daily time-series of exposure for each site.

Model building: predicting presence/absence of giant kelp
The data collected in the fine-scale analysis were used to construct binomial models to predict the presence of giant kelp canopy cover across the state as a function of SST, the state of ENSO and exposure conditions. Given the nature of the data, the model was designed to assess whether there was any overall consistency in the presence or absence of giant kelp forest habitat at regional scales associated with these physical variables; importantly, it was not appropriate with these data to attempt a model to correlate individual spikes or drops in canopy extent to particular disturbance events or temperature ano malies. Several modelling techniques were investigated, including generalised linear models (GLMs), least absolute shrinkage and selection operator (LASSO; Tibshirani 1996) and random forest (Brei man 2001). Each of these methods produced very similar results, but the LASSO was considered the most appropriate because it is intended for use with large numbers of correlated variables and produces much simpler and more stable results compared to other variable selection techniques (Tibshirani 1996). In the presence of strongly correlated predictors, the coefficients of a GLM are only weakly determined. To address this, the LASSO nor-malises the predictors so they are of comparable scale and adds a penalty of the form: (2) to the usual GLM negative log likelihood, where β i are the regression coefficients for the rescaled predictors and λ is a regularization parameter. This has the effect of shrinking the fitted coefficients towards zero at the expense of model fit; the larger the regularization parameter the greater the shrinkage, with weakly determined coefficients more strongly impacted. In this way, LASSO produces sparse fits, and has been shown to be robust in the presence of correlated predictors (Oyeyemi et al. 2015).
Random forests are also robust in the presence of correlated predictors. A random forest draws many bootstrap samples from the original data, and then fits a classification tree to a random subset of the predictors for each bootstrap sample (Friedman et al. 2001). Predictions for the model are made by averaging the predictions from the many bootstrap samples. Random forests are robust to the presence of correlated predictors because the correlated predictors are spread at random across the many trees so that each predictor contributes to the final prediction.
For these reasons, while all 3 models were fitted, we only report on the results of the LASSO and random forest.

Model parameters
Model parameters included daily mean SST estimates for each site (as in Section 2.3.1), monthly state-wide SOI (as in Section 2.3.2) and a 'storm' index. For exposure, and with the expectation that only the more extreme energy events would impact giant kelp populations (Jones et al. 2015), the WEI time-series was used to develop a storm index. The index was generated by determining how much the wave energy at a particular site exceeded a certain threshold, which was defined as a 'storm' event. The threshold was determined as the 90 th percentile of the WEI time-series for that site. If the daily wave energy was below the threshold, it was given a value of zero (i.e. not a storm); thus storms were described by their intensity (the magnitude of the WEI above the 90 th percentile threshold) and duration (number of contiguous days above the threshold). The reasoning behind using a site-specific threshold is that giant kelp morphology is plastic and to some extent correlated with site-specific conditions (Demes et al.

Model construction
The presence of a surface canopy of giant kelp on any given day is reliant to some extent on environmental conditions prior to that day. Therefore, statistical modelling used SST, exposure and SOI at several time lags. For each observation of kelp canopy, the average SST and SOI were calculated over contiguous 90 d blocks preceding the date on which the image was taken. Mean values were used because our questions relate to long-term trends in average conditions, as opposed to specific acute events. For the exposure index, values were summed to give a number that represented both the severity and duration of storm activity in each contiguous 90 d block. Previous research has shown that the importance of ocean conditions in affecting giant kelp cover and biomass extends at least as far back as 3 yr prior to observation , and thus these lagged terms were calculated for the 4 yr prior to each satellite observation being captured. For example, if kelp canopy was observed on 1 September 2000, the value for the first SST lag (3 mo) was the average of the daily SST from 1 June to 1 September 2000, the second SST lag was the average from 1 March to 1 June 2000 (i.e. also a 3 mo period, but beginning 6 mo prior to the day of observation of kelp surface canopy cover), and so on. The number of lags totalled 15 for both SST and exposure variables (= 4 yr), while historic data for SOI allowed for 20 lags (= 5.25 yr).
The potential for multiple correlated predictors in this model is high, and thus the best we may expect from our model is that it assigns larger coefficients to variables that are correlated with variables that are truly predictive of the outcome. While this may be the case, our view is that the model is robust and useful if it includes relevant predictors, reproduces observed responses well, can be explained by a set of physical processes, and where the interpretation is consistent with other evidence and observations in the current literature.
Due to a large proportion of zeros (where an observation was possible but no canopy was observed), canopy observations were converted to a binary presence/ absence classification, and the LASSO analysis took the form of a logistic regression. All variables were normalised prior to fitting the model sequence using the inbuilt function of the 'glmnet' package in R. Coefficients are reported on the original scale. LASSO solutions were computed for a range of values of λ, and the best-fit and the simplest models were selected using cross-validation by year. The best-fit model (λ.min) refers to the model where the value of λ gives the lowest cross-validation error, while the simplest model (λ.1se) refers to the model that has the fewest predictors but where the value of λ gives a cross-validation error within 1 SE of the best-fit model. The simplest model thus cannot be distinguished from the best-fit model in terms of error, given the uncertainty associated with the cross-validated estimate of error of the best model. Random forest modelling used the same set of lagged predictors and binary response of kelp surface canopy. Three thousand trees were generated. Variable importance plots were used to assess the importance of each predictor variable, and partial dependence plots were used to provide graphical representation of the marginal effect of the 16 most important variables on the probability of canopy absence.
Cross-validation by year was chosen because one of the issues inherent in collating this type of data is that observations are not wholly independent, and a degree of temporal auto-correlation exists. Once a giant kelp canopy has established, it is likely that it will still exist in observations over the following months, and in many cases, in subsequent years. Cross-validating by year enabled verification that the models were robust to auto-correlation.

State-wide extent of canopy cover
The state-wide time-series shows that the distribution of dense surface canopies varies spatially, with the most extensive canopies occurring in the southeast region of the state. A state-wide peak in canopy cover was driven predominately by dense surface canopies occurring in the Bruny bioregion (Fig. 2B), which accounted for 274.6 ha (65%) of the maximum cover observed through the time-series. Canopy ex -tent in the Davey bioregion reached a maximum of 70.8 ha in 1997−1999, while the Freycinet and Franklin bioregions contributed relatively little to the statewide total, with a maximum of < 45 ha being recorded in each over the entire time-series. The Boags and Otway bioregions had < 5 ha collectively at any one time, while no cover was observed in the Flinders bioregion; due to this low or absent cover, the Boags, Otway and Flinders bioregions were excluded from further analysis.
This time-series also showed that the total extent of giant kelp canopies across Tasmania varied significantly through time ( Fig. 2A). The first 2 decades within the time-series were characterised by periods of higher canopy cover, while over the last decade there was relatively little surface canopy of giant kelp. Peaks in the total state-wide extent of surface canopy were observed over the time periods from 1994− 1996 and 1997−1999, with a maximum area of 422.2 ha reached during 1997−1999. The total state-wide extent during all other periods was typically <100 ha. By 2015, total cover had fallen to < 10 ha.

Trends in abundance of canopy cover
The broad patterns from the fine-scale time-series re flected those obtained in the state-wide analysis (Fig. 3). Large spatial and temporal variation in cano py cover is clear, and the most extensive canopy areas were observed in the Bruny bioregion, followed by Davey, Freycinet and Franklin. With few exceptions, giant kelp canopies were most abundant and extensive during the first 2 decades of the time-series, with a period of consistently high canopy cover occurring between 1995 and 1999. A state-wide shift in canopy extent is apparent, with much less canopy observed after the year 2000. This is particularly evident in the Freycinet and Bruny bioregions, where a maximum of only 17 ha total was ob served in any given year over the latter 16 yr of the time-series.
In general, short periods of surface canopy growth and expansion are followed by rapid decrease or disappearance in canopy area occurring over monthly and annual time-scales. Major peaks in surface cano py extent are often simultaneous between sites and bioregions, and occurred in 1987, 1995, 1998 and 1999, with lesser peaks observed in 1993, 2005 and 2006. 3.1.3. Canopy cover estimates are not sensitive to tide height A 1-sample t-test on paired images demonstrated that the percentage change in canopy area did not differ significantly from zero when there was a change in tide height category (t 33 = −0.51, p = 0.63). A 2-sample t-test showed that the average percentage change in canopy area did not differ significantly between pairs where there was a difference in tide height category and pairs where there was no difference in tide height category (t 87 = 1.13, p = 0.26).

SST
All bioregions showed distinct annual cycles of SST, which were superimposed upon longer-term cycles coinciding with those of the ENSO (Fig. 4).

ENSO
The oscillations of the ENSO climate cycle (i.e. time to cycle through an El Niño and a La Niña), ranged from 2 to 6 yr, and both positive (La Niña) and negative (El Niño) extremes were observed over the time-series.   Five El Niño events occurred between 1986 and 2016, with particularly strong occurrences in 1987− 1988, 1991−1992−1998. Three consecutive events occurred from 1991− 1995, and individual events occurred in 2002. La Niña events oc curred in 1988− 1989, 1998− 2001− 2012. The events in 1988− 1989−2012 were the strongest.

Exposure
Different bioregions demonstrated different average wave exposures, with the Bruny bioregion being the most protected and sites in the Franklin bio region being the most exposed (Fig. 5). The most exposed bioregions (Franklin and Davey) had the most pronounced seasonal signal in exposure, although no long-term trends were clear (Fig. 5). Variation in exposure among sites within bioregions was also high, and some sites were clearly more exposed than others to the same high wave energy events (see Text S1 and Figs. S1−S4).

Modelling canopy presence using environmental correlates of canopy change
Binomial models were constructed using LASSO with the purpose of predicting the presence or absence of giant kelp canopy across the state based on historic SST, exposure and SOI conditions. SST variables dominated both the best-fit and simplest LASSO solutions. In both models, SST variables had coefficients 1−2 orders of magnitude greater than SOI variables, and 2−4 orders of magnitude larger than exposure variables (Table 1). In the best-fit model, the 3 most important variables were SST at lags of 15, 36 and 39 mo, with coefficients of −0.554, −0.400 and −0.154, respectively ( Table 1). The next highest coefficients were −0.059, −0.053 and −0.046 for SST lags of 27, 48 and 45 mo, respectively.
In the simplest model, SST lags of 15 and 36 mo, and SOI lag of 63 mo, were the only variables remaining in the model. SST 15 and 36 had coefficients of −0.215 and −0.128 respectively, and SOI 63 had a much smaller value of −0.009.
A boxplot of predicted vs. fitted values shows that these models performed well (Fig. 6). For a site where canopy was absent, in 75% of cases the best fit model had a probability of predicting presence of 43%, while for a site where canopy was present, the best fit model had a probability of predicting presence of 80%.

DISCUSSION
Population dynamics in giant kelp are influenced by a complex array of abiotic and biotic variables. Reproduction, recruitment, growth and mortality are all influenced by water nitrate concentrations, temperature, wave disturbance, light levels and grazing (Graham et al. 2007), as well as by human impacts such as harvesting (Vásquez 2008) and waste pollution (Foster & Schiel 2010). The dependence of growth and mortality on these environmental conditions results in a species that demonstrates a high propensity for both rapid decline and rapid recovery. The forests of Tasmania are no exception, and the results presented in both our Observed class Fitted probability of presence λ.1se λ.min Fig. 6. Distribution of fitted probabilities from the LASSO regression within each observed class (present or absent) for the best-fit model (λ.min), in which the shrinkage parameter (λ) gives the minimum average cross-validation error, and for the simplest model (λ.1se), in which λ gives a cross-validation error within 1 SE of the cross-validated error of the best-fit model. The lower and upper limits of the box represent the 25 th and 75 th percentiles, respectively (the interquartile range, IQR), while the lower and upper whiskers extend no further than the value that is closest or equal to 1.5× the IQR below/ above the box limits. Outliers are shown as circles, and the centre bold line represents the median (50 th percentile)  Table 1. Coefficients of variables in the best-fit model (λ.min), in which the shrinkage parameter (λ) gives the minimum average cross-validation error. Cells marked with a dash (−) indicate those variables whose coefficient was shrunk to zero and therefore are not considered to be a predictor of Macrocystis pyrifera canopy presence or absence fine-scale and state-wide analyses show a high degree of spatial and temporal fluctuation in cover. Notably, this study records a clear shift through time in the abundance and distribution of giant kelp in Tasmania, and when taken in the context of previous work (Sanderson 1987, Seacare 1999, 2019, Johnson et al. 2011, Steneck & Johnson 2014, an overall and persistent state-wide reduction in surface canopies is clear, with the data in the present work representing the tailend of the decline. Our models suggest a strong association of kelp forest disappearance with SST. These results are in contrast to other studies of Macrocystis pyrifera at regional scales. Reed et al. (2011) established that net primary production (as a function of biomass) in giant kelp forests was determined more by regional differences in wave disturbances than by nutrient supply (correlated with SST) and/ or grazing, while Cavanaugh et al. (2011) found no net change in biomass of giant kelp in the Santa Barbara Channel over the 25 yr period from 1984. Similarly, Friedlander et al. (2020 found no long-term trends in kelp canopy area from 1998− 2020 at the southernmost tip of the South American continent. These differences highlight the dynamic nature of giant kelp populations at a global scale, and contribute to understanding how global, climate-driven changes can interact with regional influences to alter the structure and function of coastal ecosystems (Krumhansl et al. 2016).

Importance of SST
The importance of SST in driving population dynamics of giant kelp is well established (reviewed by Schiel & Foster 2015). It is difficult, however, to unequivocally separate the effects of SST from the numerous factors with which temperature is correlated in natural systems, e.g. nutrient loading (Zimmerman & Kremer 1984) and, in some cases, salinity (Ridgway 2007a). Arguably the most important of these relationships in the context of giant kelp is the common inverse relationship between temperature and nitrate (Zimmer man & Kremer 1984). M. pyrifera sporo phytes have a limited capacity to store nitrogen (Gerard 1982), a nutrient essential for photosynthesis, and reduced nitrate levels (<1 μM ambient nitrate) can severely impact growth (Brown et al. 1997), reproduction (Reed et al. 1996), development of early life-history stages (Carney & Edwards 2010) and recruitment success (Ladah & Zertuche-González 2007) in this species.
Knowledge of the independent effects of temperature and nitrate availability on M. pyrifera comes from manipulative experiments. North & Zimmerman (1984) investigated the effect of artificial nitrogen fertilisation on mature M. pyrifera exposed to a range of temperatures above its typical range (>18°C) and found that sporophytes lost their canopy at temperatures of 18−23°C unless nutrients (nitrate and phosphate) were added, in which case the canopy thrived. Similarly, Hernández-Carmona et al. (2001) found increased recruitment and survival of juveniles transplanted to nutrient-fertilised areas relative to unfertilised controls during periods of high temperatures (up to 28°C). Conversely, recent laboratory experiments with juvenile M. pyrifera sporophytes in Tasmania using a fully factorial design indicated greater sensitivity to warming (within ob served ranges) than to nutrient depletion (Mabin et al. 2019). Certainly higher temperatures can lead to higher metabolic rates (Burdett et al. 2019), in creas ing the demand for and use of nitrate supplies (Staehr & Wernberg 2009, Buschmann et al. 2014. The overall picture to emerge from these observations is that in worst cases, recruitment, growth and survival of M. pyrifera sporophytes are negatively impacted by warm water outside the typical range experienced by a given population, but that the deleterious effects can be ameliorated to some extent by elevated nutrient levels. This is consistent with current observations in eastern Tasmania, where nitrate concentrations in surface waters are often <1 μM, and typically reach 0 μM (i.e. undetectable) in late summer (Rochford 1984), and remnant M. pyrifera individuals are typically only associated with a local source of nutrient input, e.g. sewage outfalls and stream or river mouths.

Broad-scale variability
The effect of SST around Tasmania appears to be expressed predominately through ENSO cycles, and specifically El Niño events. Our data showed that all peaks in canopy cover corresponded with El Niño events, which occurred in 1987−1988, 1991−1992, 1993− 1994, 1994−1995, 1997−1998, 2002−2003, 2006− 2007, 2009−2010 and 2015−2016. In Australian marine environments, an El Niño event is characterised by cooler than average water temperatures and reduced storm activity (Harris et al. 1988). Weakened east− west trade winds reduce the supply and transport of warm, nutrient-poor surface waters southwards along Australia's east coast to Tasmania, and allow for increased influence of colder (and more nutrient-rich) Southern Ocean water masses (Harris et al. 1988). Together these create a favourable envi-ronment for the growth and development of giant kelp. Larger than normal decreases in SST occurred during the El Niño events of 1987−1988, 1991−1992and particularly that of 1997−1998.au/climate/ enso/ enlist/), and this is reflected in some of the most extensive increases in canopy area.
Inter-annual variability of giant kelp populations is closely linked to the ENSO cycle around the globe , Friedlander et al. 2020. The most commonly documented impacts are widespread mortalities and canopy loss in the eastern Pacific that result from increased wave activity and temperature regimes associated with El Niño conditions, i.e. the opposite effect to that which occurs in Australia during El Niño events (Zimmerman & Robertson 1985, Hernández-Carmona et al. 2001, Edwards & Hernández-Carmona 2005. The period from 1993− 1999 in Tasmania saw 5 years with cooler, nutrientrich El Niño conditions, which provided a consistent and favourable environment for growth, canopy development, reproduction and recruitment, and resulted in an extended period of dense and relatively abundant canopies around the state (al though note that in eastern Tasmania this represented only ~30% of the peak canopy cover since 1946; see Steneck & Johnson 2014). Outside of El Niño periods, however, dense surface canopies have become increasingly sparse, as background water conditions become warmer and more nutrient-poor.

Regional variability
At a bioregional level, the effect of SST on the cover of dense surface canopies is consistent with ongoing changes to coastal waters that are occurring in the region. The east coast of Tasmania is a hotspot for oceanographic change on decadal time scales, with an estimated rate of ocean warming at 3−4 times the global average (Ridgway 2007a). This is due to the increased influence of EAC water in eastern Tasmania (Oliver & Holbrook 2014), and reflected in marked shifts in the temperature and salinity regimes (Ridgway 2007a). The EAC is a warm, relatively saline (Ridgway 2007a) and nutrient-depleted current with nitrate levels <1 μM (Harris et al. 1987), which is limiting to M. pyrifera growth.
This study represents the tail-end of a long-term and persistent decline in surface canopy area, particularly on the east coast of the state. In this context, the data in the present work reveal a near-complete loss of canopy cover in the Freycinet and Bruny bioregions. As the northern-most bioregions on the east coast, Freycinet and Bruny are those most influenced by the EAC. Statistical analysis of SST revealed both bioregions to be warming at a rate of 0.029°C yr −1 , with temperatures in Freycinet bioregion on average 1°C warmer than in the Bruny bioregion. This long-term increase in temperature and simultaneous decreases in nitrate concentrations will have affected the growth and survival of M. pyrifera, particularly in the Freycinet bioregion, leading to poor recruitment and recovery following summer canopy losses and disturbance events. Research by Hollar smith et al. (2020) also suggests that the increase in temperature may be influencing production of eggs or progression to the diploid life-history stage in microscopic gametophytes.
Although the southern extension of EAC water is likely the cause for the long-term decline in canopy area on the east coast, it is not possible to isolate the relative importance of the different environmental factors associated with this change. Alongside the ef fects of increased temperature and decreased nutrients, there is also increased salinity. Although Busch mann et al. (2004) suggested that giant kelp popu lations may have differentiated responses, and in some cases broad tolerance, to variable salinity levels, the majority of literature has investigated this only for salinities lower than typical oceanic concentrations (< 34 ‰; e.g. North et al. 1986, Peteiro & Sánchez 2012. While it has been suggested that the salinity changes associated with the EAC are unlikely to be sufficiently large to significantly influence the population dynamics of M. pyrifera in Tasmania (Johnson et al. 2011), this may warrant further investigation.
The transport of warmer waters associated with the EAC has also enabled the establishment of a number of new species in Tasmania, including the longspined sea urchin Centrostephanus rogersii (Johnson et al. 2005, Ling et al. 2009). C. rogersii has the ability to catastrophically over-graze kelp, and as populations have increased, they have reduced large tracts of once healthy shallow rocky reef habitat to socalled 'urchin barrens' (Johnson et al. 2005). While the long-term declines of giant kelp on the east coast of Tasmania have not been directly linked to the establishment of these urchins, it is likely they have contributed to canopy losses, perhaps by preventing recovery of populations through grazing of gametophytes, recruits and juvenile sporophytes.
The detrimental impact of the EAC on giant kelp populations is also highlighted by the less pronounced shifts in the abundance and extent of sur-face canopy cover observed on the south and west coasts. Statistical analyses of SST data showed that the Davey and Franklin bioregions had average temperatures and warming trends that statistically differed from all other bioregions. In these regions, the EAC has little influence, and the oceanographic environment is instead dominated by the ZC (Franklin bioregion) and the influence of Southern Ocean water masses (Ridgway 2007b).

Importance of time-lagged effects
The time lags selected by the LASSO models as the best predictors of canopy presence/absence were SST lagged at 15, 36 and 39 mo (best-fit model; the same variables were also selected using a random forest model) and SST lagged at 15 and 36 mo and SOI at 63 mo (simplest model). Together these correspond to time lags of ~1, 3 and 5 yr.
The annual signal is presumably driven by canopy development associated with the spring growth maximum, with subsequent deterioration over summer. As elsewhere, growth of M. pyrifera in Tasmania is highly seasonal. The only thorough assessment in Tas manian waters was conducted by Cribb (1954), who observed that maximum growth rates occurred during May−September, with the densest canopies developing in September once winter storms had subsided and light levels began to rise while nitrate was relatively available. Similar results have been re corded in New Zealand, where growth rates of M. pyrifera show a significant decline in summer resulting from limiting nitrate conditions (Brown et al. 1997), and also in southern Australia, where a biomass accumulation in laminarian species has been shown to peak in spring (Fairhead & Cheshire 2004). A dense canopy that forms in winter/spring may significantly deteriorate over the latter part of summer, but given favourable conditions will rapidly re-form in the next winter/spring from subsurface sporophytes (although we note the case of annual populations in protected areas of Chile, which show a different dynamic; Buschmann et al. 2006).
The importance of SST at a 3 yr time lag is consistent with results found by Cavanaugh et al. (2011) in a similar study applied to kelp forest biomass in the Santa Barbra Channel. It also agrees with observations of post-disturbance recovery periods of ~2−3 yr in forests along the entire east Pacific coast (Graham et al. 1997, Edwards & Hernández-Carmona 2005, Edwards & Estes 2006). This lag possibly reflects the period required for recruitment and recovery following disturbance events and the mortality of entire sporophytes.
Previous research has demonstrated that conditions immediately following disturbance events play a critical role in determining community structure and dynamics through their influence on recruitment and growth of juvenile sporophytes. Tegner et al. (1997) compared growth and survival of 2 different cohorts of giant kelp in contrasting oceanographic conditions, and found that cooler and nutrient-rich conditions led to competitive dominance by M. pyrifera, while in warmer and nutrient-poor conditions, the cohort was characterised by poor growth, survival and canopy formation. These effects on community structure lasted the entire life-span of the cohort, even after oceanic conditions returned to the 'neutral' state. Thus, a disturbance event followed by favourable conditions will allow recruitment of juvenile sporophytes within a year, and over the next 2− 3 yr, for these to grow and mature to form dense canopies.

CONCLUSIONS
Our results are an example of region-specific signals of global change acting synergistically with local stressors to result in trends that may not be representative of other parts of the world (Krumhansl et al. 2016). In an extensive review of the biology and ecology of kelp forests, Schiel & Foster (2015) showed that giant kelp still occupies much of its post-glacial distribution. The long-term decline in the surface cover of giant kelp in Tasmania beginning in the 1970s and 1980s is an exception, driven largely by the extraordinary changes in ocean climate associated with changes in the behaviour of the EAC (Ridg way 2007a, Oliver & Holbrook 2014), and, more recently, exacerbated by the establishment and build-up of an introduced urchin population ( Johnson et al. 2005). Importantly, this scenario has played out at the southern limit of the range of M. pyrifera in Australia, and there are no shallow reefs further to the south within dispersal distance to act as a refuge for the species in this region. In this sense, southeast Australia in general, and Tasmania in particular, represents a geographic trap for M. pyrifera (and other cool-temperate species; Edgar et al. 1991). Given the significant decline in cover of M. pyrifera in Tasmania, which is the Australian stronghold of the species, and the problem of the geographic trap, the Federal listing in August 2012 of giant kelp as an endangered marine community type in Australia is justified.
While it may be comforting to know that the status of giant kelp is stable at a global scale (Schiel & Foster 2015), our research highlights the potential for future declines in local areas or regions at risk of ocean warming. Anthropogenic impacts are increasing in areas where M. pyrifera habitat is proximal to urban areas, so in many areas, the response of this species to the complexity of multiple stressors is difficult to predict with certainty.