DOI: 10.5061/dryad.xd2547dn1
title: Climate-change-driven cooling can kill marine megafauna at their distributional limits

description: The impacts on marine species from secular warming and heatwaves are well
 demonstrated, however the impacts of extreme cold events are poorly
 understood. Here, we link the death of organisms from 81 species to an
 intense cold upwelling event in the Agulhas Current, and show trends of
 increasing frequency and intensification of upwelling in the Agulhas
 Current and East Australian Current. Using electronic tagging, we
 illustrate potential impacts of upwelling events on the movement behaviour
 of bull sharks, Carcharhinus leucas including alterations of migratory
 patterns and maintenance of shallower dive profiles when transiting
 through upwelling cells. Increasing upwelling could result in
 “bait-and-switch” situations, where climate change expands subtropical
 species’ distribution, while simultaneously exposing climate-migrants to
 increased risk of cold-stun/mortality events at poleward distributional
 limits. This shows the potential impacts of increased cold events, an
 understudied aspect of climate-change research, and highlights the
 complexities of climate change effects on marine ecosystems.
description: <strong>Study areas</strong> Both
 study areas are located along Western Boundary Currents (WBCs): the
 Agulhas Current in southern Africa and the East Australian Current in
 eastern Australia. In southern Africa, the area spans from the
 warm-temperate Breede River estuary on the south coast of South Africa
 (34.5°S) to northern Mozambique (16°S) (Figs. 1, 5a main MS). In eastern
 Australia, the area extends from the warm-temperate Sydney Harbour estuary
 (33.5˚S) to tropical Queensland (16.5˚S) (Fig. 5b main MS). On both
 continents, the study areas were divided into two distinctive zones (Fig.
 5a, b main MS). First, an upwelling zone, comprising multiple coastal
 upwelling cells inshore of each WBC, extending from the southern limit
 equatorward to 32.4°S along the Agulhas Current and 28.6°S along the East
 Australian Current. This zone includes the Breede River and Sydney Harbour
 estuaries in South Africa and Australia, respectively (Fig. 5a, b main
 MS). The second zone, in each case, was the subtropical/tropical section
 of coastline equatorward of the upwelling zones. Both
 the Agulhas and East Australian Current closely follow the continental
 shelf, with their warm cores usually situated between 50-100 km offshore
 from the shelf edge16, 56. Note that although upwelling can occur in the
 subtropical/tropical areas to the north, it is most prominent and intense
 in both upwelling zones16, 30. In the case of the
 Agulhas Current, coastal upwelling can be caused by current forcing as
 well as by meanders in the current<sup>16, 19, 31</sup>. Such
 meanders can be associated with cold-core inshore eddies that travel
 south-ward along the shelf, often called “Natal Pulses”<sup>
 17</sup>. Current-driven coastal upwelling is most common along the
 south-east coast of South Africa, due the proximity of the shelf edge to
 the coast<sup>16, 19, 31</sup>. Upwelling-favourable easterly
 winds can also contribute to coastal upwelling
 here<sup>16</sup>. As the shelf widens south of Port
 Elizabeth, this becomes the main mechanism of coastal upwelling, with a
 weakening influence of the current. Similarly, along the East Australia
 Current, current and eddy-driven upwelling mostly occurs along the narrow
 shelf from 24˚S to 32.5˚S, after which the current separates and
 wind-driven upwelling becomes the most prominent mechanism<sup>31,
 32, 56</sup>. <strong>Collection of
 environmental variables </strong> On March 2nd
 2021, dead marine organisms began washing up on beaches inshore of the
 Agulhas Current,  between Port Elizabeth to ca. 100 km north of East
 London (Fig. 1a main MS). Specimens were collected and cataloged by Greg
 Brett in association with the East London Museum (Supplemental Table 2).
 To determine whether rapid-temperature drops associated with upwelling
 could be responsible for the deaths and, if so, to qualify characteristics
 of lethal upwelling events, we investigated sea-surface temperature (SST)
 from satellite imagery, in-situ temperatures at 30 m depth and wind speed
 and direction. We used the Aqua Modis sea-surface
 temperature data set from the Moderate Resolution Spectroradiometer on the
 AQUA NASA satellite, a daily, gridded level-3 data set at 0.04166°
 resolution, with data availability from 2003-2022. To analyse the duration
 as well as maximum spatial extent and rate of temperature declines for
 potential upwelling before, during and after the dead specimens washed up
 on local beaches, we generated mapped images in ArcGIS based on
 NetCDF-files downloaded from the OceanColor webserver (https://oceancolor.gsfc.nasa.gov/l3/). Satellite-derived temperature data was complemented with in-situ temperatures collected by a temperature logger deployed 3 km offshore from Port Alfred (Fig. 1a main MS) at 30 m,, which recorded temperature at hourly intervals. Finally, wind data collected by the South African Weather Service at Port Alfred was used to determine whether upwelling-favourable easterly winds (45°-135°) occurred before and during the period dead marine organisms washed up on local beaches. Wind speed in m/s and wind direction in degrees was collected at 6-hour intervals daily between 8 am and 8 pm local time. <strong>Modelling trends in upwelling characteristics over time</strong> To investigate whether there has been an increasing trend in the proportion of upwelling-favourable winds on the south and south-east coast of South Africa, we used wind data collected by the South African Weather Service from three stations within the upwelling zone. Note that a similar analysis has previously been done for eastern Australia<sup>35</sup>. Daily wind speed in m/s and wind direction in degrees was collected at 6-hour intervals between 8 am and 8 pm local time. Data was available from November 2001 until March 2022 at Port Alfred, from October 1992 until March 2022 at Port Elizabeth, and from October 1988 until March 2022 at Plettenberg Bay. Based on the orientation of the coastline in the upwelling zone and wind direction in degrees, the data set was divided into upwelling-favourable winds (north-east to south-east (45°-135°) and non-upwelling favourable winds (all other directions)<sup>57</sup>. To estimate the daily prevailing wind speed and direction, data was averaged per day. Since we aimed to investigate potential impacts of increasing trends in upwelling on mobile marine species at their distributional limit, only wind data between October and April was included in the models. Based on bull shark tracking data from this study, this was the period during which animals tagged in the Breede River and Sydney Harbour estuaries were mainly present in the upwelling zone. Furthermore, this period also encompasses the austral summer months which are the main period of upwelling activity inshore of the Agulhas and East Australian Currents<sup>16, 19, 30, 34</sup>. Thus, annual data from October to April were combined to represent corresponding upwelling seasons. As a result, we had data for 21 upwelling seasons at Port Alfred, 30 at Port Elizabeth and 34 at Plettenberg Bay. Although timeseries at Port Alfred were shorter than 30 years, which is the recommended length for trend analysis<sup>21</sup>, we included this site to inspect a potential trend and compare it with Port Elizabeth and Plettenberg Bay sites. We then calculated the proportion of days during which upwelling-favourable winds occurred in each upwelling season. The resulting values were modelled against time (upwelling season) using generalised linear models (GLM) with binomial error distribution. Model fit and diagnostics, including tests for temporal autocorrelation were evaluated using the Dharma package in R. Only trends in the South African study areas were modelled as those in the east Australian coastline has previously been examined<sup>35</sup>. In order to investigate whether the number, mean intensity (in °C) and the mean rate of onset (in °C/day) of potential upwelling events has increased in upwelling cells inshore of both WBCs, we used the long-term, interpolated, reanalysed, daily Global Ocean OSTIA Sea Surface Temperature data set<sup>58</sup>. This gridded, level-4 data set covers the period from 10/01/1981 until present at a resolution of 0.05°, combining satellite data as well as in-situ data. We chose this data set because of long-term data availability of 41 years, the high standard of quality control and overall accuracy and applicability<sup>59</sup>. Although level-4 data sets can underestimate thermal signatures of upwelling zones due to warm biases, the OSTIA data set has been shown to perform better than other long-term data sets<sup>36</sup>. Furthermore, this bias is more pronounced in Eastern Boundary Current compared to Western Boundary Current upwelling systems<sup>33, 59</sup>. Next, we extracted sea-surface temperature data for prominent upwelling cells inshore of each WBC. This was off Port Alfred, Port Elizabeth and Plettenberg Bay for the Agulhas Current and Port Stephens, Diamond Head and Ballina for the East Australian Current<sup>16, 19, 30, 34, 37, 47</sup>. We then divided the shelf into inshore (0-15 km from shore) and midshelf (15-30 km from shore) sections. For each station, a pixel was randomly selected from inshore and midshelf sections from the gridded OSTIA SST data set along a transect (Fig. 7 main MS). Each transect extended from the coastline, across the shelf in a direction perpendicular to the shelf edge. Choosing inshore and midshelf pixel allowed us to compare trends across areas of varying shelf width where each WBC would influence upwelling differentially. Additionally, as a control, we selected an offshore pixel within the warm core of each current along the same transects for each station. No upwelling occurs within the offshore core of the currents and possible cold events would be due to temperature anomalies other than upwelling, hence no trends should be evident. We used the R package heatwaveR to detect cold events in each randomly selected pixel over time<sup>60</sup>. The algorithm implemented in this package is based on the hierarchical definition of extreme temperature events in21 and follows a similar methodology to characterising such events. In short, marine extreme temperature events are defined as discrete, anomalously warm or cold events with characteristics such as intensity, rate of onset, and spatial extent<sup>21</sup>. For extreme cold events, these characteristics are calculated in cases when temperatures fall below the 10th percentile threshold, which is based on a baseline climatology at least 30 years in duration to avoid impacts of decadal variability. Thus, data sets with long-term data availability, such as the OSTIA SST data set, available from 01/01/1981-present, are preferred<sup>21</sup>. Before calculating the baseline climatology at each pixel, we detrended the daily OSTIA sea surface temperature data as long-term trends in the climatology can affect detection of extreme temperature events using the detrend() function61. Within the detect_event() function of the heatwaveR package, we set the minimum duration for cold events to one day, in order to include short-term upwelling events which are common along both WBCs57. From the output of this function, we summed all discrete cold events each year from October to April, the time period tagged sharks were mainly present in both upwelling zones. We also calculated the annual mean intensity (in °C) and mean rate of onset (°C/day) for events from October to April. Visual inspection of the output suggested that trends in either variable remained similar using detrended and non-detrended data. We used generalised linear models to investigate trends in the number, mean intensity and mean rate of onset of upwelling events at each pixel between 1981 and 2022, with upwelling season (aggregated annually between October and April) as the independent variable. The number of upwelling events was modelled assuming a Poisson error distribution with a log-link function, while mean intensity and mean rate of onset were modelled assuming a Gaussian error distribution. Model fit and diagnostics were inspected using the Dharma package. When significant non-linearity and/or quantile deviations were detected, we fitted generalised additive models (GAMS) (Supplemental Tables 6 and 7), which were also inspected using Dharma and the gam.check() function in the mgcv package. Finally, as decadal variability and biases associated with L4-satellite SST products can confound trends in upwelling signals, we utilised temperature logger data from four locations within the upwelling cells of the Agulhas current (deployments from 1991 to 2022) and three different satellite SST products to validate our analyses above. Loggers were Star OddiTM temperature loggers with accuracy of ± 0.05°C located at Ystervarkpunt (78 km east of the Breede River), Mossel Bay (127 km east of the Breede River), Tsitsikamma (37 km east of Plettenberg Bay) and Mosterts Hoek (85 km west of Port Elizabeth) and serviced every 1-2 years (for deployment details see Supplemental Table 5). Satellite SST products used were the L4 OSTIA SST dataset from our original analyses (5 km resolution), the L4 CCI SST product (5 km resolution) and the L3 Pathfinder AVHRR night and day products combined into a single product (5 km resolution). Warming/cooling trends in SST products and logger data were estimated using the Theil-Sen estimator. The estimator is nonparametric and hence more robust than the parametric least-squares regression as it is resistant to outliers and does not draw from any probability distribution. Theil-Sen uses the median to estimate the slope as opposed to least-squares which makes use of a weighted mean for slope estimation. Spatial SST trends were calculated using each grid point within the designated area of the gridded SST product. Additionally, the SST and logger data trends were calculated using the anomalies to avoid the complicating effects of seasonality on trend values. Trend analyses for all logger locations and SST products ranged from June 1991 to August 2022. <strong>Animal tagging and tracking </strong> Tagged bull sharks were part of previously published data sets from southern Africa and eastern Australia and were acoustically tracked in coastal receiver arrays. For details about tagging methods and receiver arrays see27, 39, 62 , Supplemental Table 8). In short, in Australia, 25 bull sharks were internally tagged in Sydney Harbour with InnovaSeaTM V16TP tags that also recorded depth and temperature at the time of detection (tag nominal delay 30-90 s). Out of 41 bull sharks acoustically tagged in southern Africa as part of a separate study 27, 62, 15 individuals were also equipped with temperature and depth recording pop-up satellite archival transmitters (Wildlife ComputersTM Mk10-PSAT and miniPAT tags) and included again in this study (Supplemental Table 8). Five of these sharks were tagged inside the warm-temperate Breede River Estuary, and 10 at the subtropical Pinnacle Reef. PSAT-tags archived depth, temperature, and light level data at pre-programmed intervals of three seconds (n=1), five seconds (n=11) and 25 seconds (n=3). Additionally, binned data at 24-hour intervals was collected for temperature (9-11°C, 11-13°C, 13-15°C, 15-17°C, 17-19°C, 19-21°C, 21-23°C, 23-25°C, 25-27°C, 27-29°C, 29-31°C, 31-33°C) and depth (0-5m, 5-10m, 10-20m, 20-30m, 30-40m, 40-50m, 50-100m, 100-200m, 200-300m, 300-400m) for all PSAT-tags. Tags were released from animals after pre-programmed deployment durations (90-360 days). Tags then floated to the surface and transmitted data to the Argos satellites. Out of the 15 PSAT-tags deployed, 10 were physically recovered by the authors via boat, providing full archived temperature and depth data (one from a bull shark tagged at the Breede River and 9 from Pinnacle Reef). Transmitted data was initially decoded using the manufacturer’s software (WC-DAP 3.0, Wildlife Computers, Inc., Redmond, WA). The Wildlife ComputersTM GPE3 State-Space model was used to generate estimates of most likely locations from acoustic detections and PSAT-tag light level data using a swimming speed of 0.7 m/s64 and to create movement tracks for each shark (Supplemental Figs. 19-23). Then, in the statistical software R63 each shark location was matched with depth and temperature data collected by the PSAT-tags to presence in the two regions of interest: The upwelling zone (distinguishing locations inside the Breede River and outside), and the subtropical/tropical areas past the upwelling zone (Supplemental Figs. 19-23). For example, if an acoustic detection or geolocation estimate occurred in the upwelling zone (Fig. 5 main MS), the corresponding depth and temperature value at this timestamp was assigned to the upwelling zone. Note, this was only done for sharks tagged in the Breede River, as animals tagged at Pinnacle Reef remained  in subtropical/tropical areas with limited upwelling. Similarly, for sharks acoustically tagged in Sydney Harbour, detections at receivers were assigned to either being in the upwelling zone, (distinguishing detections inside Sydney Harbour and outside), or past the upwelling zone (subtropical-tropical Queensland) (Fig. 5b main MS).  Studies from multiple geographical locations have shown that bull shark abundance declines significantly when seasonal temperatures fall below a long-term thermal affinity threshold of 19°C38, 39, 40, 41. To further test this threshold and investigate thermal preferences for bull sharks, we established a tentative thermal range, using bull shark occurrence records in southern Africa and eastern Australia from the Global Biodiversty Information Facility (GBIF)65in combination with our acoustic detections. For this, we filtered combined bull shark occurrences to daily, verifed human and machine observations and spatio-temporally matched each occurrence with a pixel from the AQUA Modis L3-satellite temperature data set. Similar analysis using occurrence data was also performed for three megafauna species that have had fatalities in severe upwelling events, the whale shark Rhincodon typus manta rays, Mobula birostris and Mobula alfredi. For South African sharks for which PSATs were retrieved and full archival data was available (Supplemental Table 15), we plotted continuous temperature and depth over time, while indicating the occupied zone. Binned data from all Breede River bull sharks (full archival and partial data) was used to calculate percent of time at depth and temperatures to compare this in both upwelling and non-upwelling zones Pinnacle Reef sharks were detected only outside of upwelling cells). The maximum duration in consecutive hours, for which the Breede sharks were experiencing temperatures below 19°C for each zone was also calculated. Similarly, we binned temperature and depth data recorded from acoustic detections of bull sharks tagged in Sydney Harbour to determine percent of time at depth/temperature. Acoustic detections at depth/temperature were also plotted over time for sharks tagged in Sydney Harbour and indicated in which zone they occurred. Lastly, data from temperature loggers deployed outside the Breede estuary at 20 m and off Port Alfred at 30 m, as well as the loggers used for Theil-Sen estimation (see previous section)  were used to spatio-temporally match temperature and depth recorded by PSAT-tags of the Breede River sharks with upwelling events at nearby loggers based on shark acoustic detections/PSAT-geolocations. This was plotted in combination with a detection timeline from combined acoustic and PSAT-tag locations to indicate when and where PSAT-tagged sharks encountered active upwelling cells. Temperature data for outside Sydney Harbour was recorded by a mooring located 14 km offshore, between 12-20 m depth, which was used to plot upwelling events with acoustic shark detections. <strong>Method References</strong> 56.          Archer MR, Roughan M, Keating SR, Schaeffer A. On the variability of the East Australian Current: Jet structure, meandering, and influence on shelf circulation. Journal of Geophysical Research: Oceans 2017, 122(11): 8464-8481. 57.         Abrahams A, Schlegel RW, Smit AJ. A novel approach to quantify metrics of upwelling intensity, frequency, and duration. PLoS One 2021, 16(7): e0254026. 58.          Good S, Fiedler E, Mao C, Martin MJ, Maycock A, Reid R, et al. The current configuration of the OSTIA system for operational production of foundation sea surface temperature and ice concentration analyses. Remote Sensing 2020, 12(4): 720. 59.          Yang C, Leonelli FE, Marullo S, Artale V, Beggs H, Nardelli BB, et al. Sea surface temperature intercomparison in the framework of the Copernicus Climate Change Service (C3S). Journal of Climate 2021, 34(13): 5257-5283. 60.          Schlegel RW, Smit AJ. heatwaveR: A central algorithm for the detection of heatwaves and cold-spells. Journal of Open Source Software 2018, 3(27): 821. 61.          Schlegel RW, Oliver EC, Hobday AJ, Smit AJ. Detecting marine heatwaves with sub-optimal data. Frontiers in Marine Science 2019, 6: 737. 62.          Daly R, Smale MJ, Cowley PD, Froneman PW. Residency patterns and migration dynamics of adult bull sharks (Carcharhinus leucas) on the east coast of southern Africa. PloS one 2014, 9(10): e109357. 63.        R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/ 64.          Brunnschweiler JM. Sharksucker–shark interaction in two carcharhinid species. Marine Ecology 2006, 27(1): 89-94. 65.         Stuart-Smith, R. D., Edgar, G. J., Barrett, N. S., Kininmonth, S. J., & Bates, A. E. (2015). Thermal biases and vulnerability to warming in the world’s marine fauna. Nature, 528(7580), 88-92.
description: # Climate-change-driven cooling can kill marine megafauna at their
 distributional limits Bull Shark PSAT Data and example code for modeling
 of upwelling variables: ## Description of the Data and file structure
 These data include two different data sets Data set 1 contains 5 csv.
 files with raw depth and temperature data collected by PSAT-tags
 externally attached to bull sharks tagged in the Breede River in South
 Africa (n=5, 1 per shark). For this study, recorded depth (in m) and
 temperature (in °C) were the main variables used. For Shark ID1, ID2, ID4
 and ID5 PSAT-tag series data was available. This includes the columns
 DeployID and PTT which are each sharks individual identifier number,
 DepthSensor which is the depth sensor resolution in m, a Source column
 which describes the source of the data, in the case of our PSAT-tags this
 is by Transmission, an Instrument column which describes the tag type, in
 our case MiniPAT, Day and Time columns provided in UTC. The Depth column
 provides Depth in m, while the DepthRange column provides the depth range
 experienced by each individual shark during each time bin (in our case 10
 min interval). Similarly, the Temperature column provides temperature
 values in °C with the TemperatureRange column indicating the range of
 temperatures are shark experienced in a binned time interval of 10 mins.
 Activity and ARange columns were not used for this study. LocationQuality,
 Latitude and Longitude are not given in Series data and will populated
 with NA. Lastly, for SharkID3 full Archival Data was available, thus the
 data structure for this individual differs slightly. Here, the Time column
 provides a combined date and time value in UTC, with Depth and Temperature
 given in meters and °C respectively. Lightlevel data is used to determine
 geolocation for an animal at a given timestamp by measuring light
 intensity in W cm^-2. Columns Ax, Ay, Az, Dry, Aux, One Minute Light
 Level, Smoothed Light Level and Events are not applicable to this study
 and can be ignored. Combined tracks for each shark tagged in the Breede
 River are given in supplemental figures 18-22 to match with
 depth/temperature values. All rows containing NA means that tags were
 unable to measure and/or transmit to a satellite the given value at that
 time stamp. Data set 2 includes example code for how we modelled number of
 upwelling events, mean intensity of upwelling events and mean rate of
 onset of upwelling events over time (see methods). For this, we also
 included a Net_CDF file for an example location which includes a column
 for sea surface temperature and Date. All other data used in the MS is
 publicly available at: Both the OSTIA and the CCI satellite SST product
 were accessed at
 [https://marine.copernicus.eu/](https://marine.copernicus.eu/). The AVHRR
 Night satellite SST product was accessed at
 [https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdPH53sstnmday_Lon0360.html](https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdPH53sstnmday_Lon0360.html), while AVHRR Day satellite SST product was accessed at [https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdPH53sstdmday_Lon0360.html](https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdPH53sstdmday_Lon0360.html). The L3 Aqua-Modis satellite SST product was accessed at [https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdMH1sstd1dayR20190SQ.html](https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdMH1sstd1dayR20190SQ.html). These data sets provide temperature in °C, a date and latitude and longitude in decimal degrees. NetCDF-files used to map the progression of the lethal upwelling event were downloaded from the OceanColor webserver ([https://oceancolor.gsfc.nasa.gov/l3/](https://oceancolor.gsfc.nasa.gov/l3/). Long-term wind data for South Africa can be requested from the South African Weather Service [https://www.weathersa.co.za/](https://www.weathersa.co.za/). Wind data for key dates in main manuscript is supplied in the supplemental material, including a column for wind speed in m/s and wind direction in degrees for a given date. Data associated with bull sharks tagged in Sydney Harbour (temperature in °C, depth in m and acoustic detections at a given receiver) is publicly available from the animal tracking facility at IMOS [https://animaltracking.aodn.org.au/](https://animaltracking.aodn.org.au/) (see supplemental material for tagging dates of sharks). Similarly, acoustic tracking data for bull sharks tagged in Southern Africa is available from the South African Acoustic Tracking Array Platform [https://www.saiab.ac.za/atap.htm](https://www.saiab.ac.za/atap.htm). Mooring temperature data for Sydney in °C is publicly available at the Australian Ocean Data Network on [https://portal.aodn.org.au/](https://portal.aodn.org.au/) (Mooring code: SYD100). In-situ temperature data collected by a logger at Port Alfred in in °C is supplied in the supplemental material.