Ocean warming projected to increase ecotourism opportunities to encounter iconic marine megafauna (manta rays and zebra sharks) off south-eastern Australia

Study area

South-eastern Australia contains the poleward-range edge for the distributions of manta rays and zebra sharks in eastern Australia (Couturier et al. 2011; Dudgeon et al. 2013) and, therefore, is the region where changes in annual temporal persistence of these species are most likely to occur. Consequently, habitat-suitability models were built using occurrence data for manta rays (Fig. 1a) and zebra sharks (Fig. 1b) from the Arafura, Coral and Tasman Seas located in the southern hemisphere (0–40°S) on the eastern side of Australia (130–180°E).

Fig. 1.

Study area, showing occurrence records (circular symbols) for (a) manta rays (Mobula alfredi) and (b) zebra sharks (Stegostoma tigrinum). Colours show summer maximum sea-surface temperatures (T_max). Symbol size indicates month of occurrence record.

Species and environmental data

Occurrence data were extracted from the Global Biodiversity Information Facility (see https://www.gbif.org, accessed 22 August 2023), the Reef Life Survey species occurrence database (see https://portal.aodn.org.au/search, accessed 22 August 2023), Atlas of Living Australia data (see bie.ala.org.au, accessed 28 August 2023), and dive-shop records (Supplementary Table S1). Data from all sources provided the location and date where a target-species presence was recorded. Data were screened to include only data from 2000 onward, with this period selected to match the period for the environmental data used in modelling and given that only very limited species data were available prior to this date.

Data were screened to exclude records for preserved specimens, records for zebra shark egg cases and records with no date or location-coordinate information. All citizen-science records were verified by examination of photos, to confirm correct species identification. Additionally, multiple records from the same location (within ~5-km radius of each other) and from the same week (within a 7-day window of each other) were combined to form a single record, as recommended for filtering biological data (Johnston et al. 2021) and appropriate for handling spatial and temporal autocorrelation structures within species-occurrence datasets (Brodie et al. 2015; Champion et al. 2021). The resulting screened presence data consisted of a mixture of machine detections from tagged animals carrying satellite and acoustic transmitters (38.2% manta rays, 22.0% zebra sharks) citizen-science records (61.0% manta rays, 76.4% zebra sharks) and records from scientific field surveys (0.8% manta rays, 1.6% zebra sharks; Table 1). To ensure that habitat suitability models were fitted to the same amount of data, and thus their results were subsequently comparable, 700 occurrence records were randomly sampled for both species from the resulting dataset and retained for model selection and evaluation.

To create a binomial response variable for modelling the environmental habitat preferences of manta rays and zebra sharks, pseudo-absence data were spatially and temporally randomised in water shallower than 200 m (i.e. the depth of the continental shelf break) throughout the study extent and period to characterise unsuitable habitat for each species. A consistent dataset containing a total of 10,000 pseudo-absences was used for each species, on the basis of Barbet-Massin et al. (2012) who recommend a large number (i.e. 10,000 or more) of randomly selected pseudo-absences for statistical modelling of species habitat preferences.

Multiple, dynamic ocean variables that are known to influence the seasonal distribution patterns of mobile marine species off eastern Australia (Hobday and Hartog 2014; Champion et al. 2021) were downloaded from the Copernicus Marine Environment Monitoring Service (CMEMS, see https://data.marine.copernicus.eu/products) and matched to manta ray and zebra shark occurrence and pseudo-absence data. These variables were sea-surface temperature (SST; 0.05° spatial resolution; product #010_011), eddy kinetic energy (EKE; 0.25° spatial resolution; product #008_047), sea-level anomaly (SLA; 0.25° spatial resolution; product #008_047) and chlorophyll-a concentration (CHL; 0.04° spatial resolution; product #009_082).

Collinearity among ocean variables was assessed using Spearman rank correlation coefficients and variance inflation factors (VIFs). Spearman coefficients for pairwise comparisons among ocean variables did not show any collinearity (i.e. values for Spearman’s r were <0.5 or >−0.5). Furthermore, VIFs associated with all ocean variables were <1.5, indicating that collinearity among these variables was unlikely to affect model quality (Zuur et al. 2007).

Habitat-suitability modelling

Generalised additive mixed-effects models (GAMMs) were used to parameterise habitat suitability preferences for manta rays and zebra sharks. GAMMs used a logistic link function to relate the binomially distributed response variable (i.e. species occurrence or pseudo-absence) to multiple dynamic ocean-state covariates. Models also included ‘calendar year’ as a random intercept term to control for intra-class correlations among data collected during the same year, which may have been temporally biased towards certain years because of unmeasured annual variation in species observational effort over the study period. To accommodate non-linear responses, penalised regression spline smoothers were fitted using generalized cross-validation as recommended, to optimise smooth functions and avoid overfitting to the data (Zuur et al. 2009). Smoothers were removed from predictor variables in favour of linear model fits if their effective degrees of freedom were approximately equal to 1, which indicates linearity with log-of-odds transformed binomial response variables. These models estimate the probability of species occurrence in space and time. However, because absolute probabilities are contingent on the ratio of occurrence to absence data (Pearce and Boyce 2006), and given that pseudo-absences were randomly generated herein, probabilistic predicted values were rescaled to an index of environmental habitat suitability ranging between 0 (unsuitable) and 1 (optimal).

A backward model-selection procedure was initially employed to identify significant predictors of manta ray and zebra shark occurrence by using all ocean variables. Following this, models underwent further refinement via a stepwise selection process guided by covariate P-values. This yielded a set of explanatory models for each species, characterised by nested combinations of covariates of decreasing complexity. Among these models, the one displaying the lowest Akaike information criterion (AIC) value was identified as the most parsimonious (i.e. the best model). A key objective of this study was to develop a time-series of monthly habitat-suitability projections for both study species under future climate change. However, future data from the latest generation of global climate models (CMIP6) produced at monthly time-steps and downscaled to a spatial resolution relevant to dive location along Australia’s eastern coast were available only for the SST model covariate (see Future projections section). Therefore, we retained temperature-only models (i.e. thermal-habitat suitability models) from the nested set of habitat-suitability models for each study species for the purpose of creating future projections of thermal suitability for manta rays and zebra sharks. Assessments of model accuracy and predictive skill (described below) were compared between temperature-only and best habitat-suitability models for both species, with the best models being associated with only minor improvements in model quality over temperature-only models.

Given that the objective of this study was to project species habitat suitability for a forthcoming period, we utilised an out-of-sample evaluation method to gauge model accuracy and predictive capability (Champion et al. 2022). This involved training best and temperature-only habitat models for individual species by using subsets of occurrence and pseudo-absence data from the timeframe spanning 2000–15 (referred to as Period 1), followed by model testing utilising data from 2016 to 2023 (referred to as Period 2). These periods were selected as they each contained ~50% of the observational data for both study species. Therefore, model accuracy was tested between two periods that have experienced historical climate-driven environmental changes, which provides an appropriate assessment of model accuracy for future periods associated with continued climate change.

A six-fold cross-validation was undertaken by randomly extracting six subsets of species-occurrence and pseudo-absence data, each containing 50 occurrence records and 1000 pseudo-absences for every species. Of these, five subsets were drawn from Period 1 for model-training purposes, whereas the sixth subset was drawn from Period 2 for model testing. Subsequently, best and temperature-only habitat-suitability models for each species were fitted to each Period 1 subset and subsequently assessed against the Period 2 subset. The outcomes of these evaluations yielded a series of confusion matrices, which describe rates of true positive or negative outcomes as well as false positive or negative outcomes, which were then utilised to quantify two metrics of model accuracy. These metrics included the mean area under the receiver-operating characteristic curve (AUC; range between 0 and 1), where values exceeding 0.8 signify good to excellent model performance (Araújo et al. 2005), and the mean true skill statistic (TSS; range between −1 and 1), where positive values indicated models surpassing chance expectations (Allouche et al. 2006). Differences in AUC and TSS values between the best and temperature-only habitat-suitability models for each species were quantified to compare model performance.

Future projections

Monthly thermal-habitat projections for manta rays and zebra sharks over the period 2015–59 were created using an ensemble of sea-surface temperature (SST) data downscaled from seven global climate models (GCMs; CMIP6; Table 2). GCMs were selected using the IPCC Working Group 1 Interactive Atlas of Regional Information interface (see https://interactive-atlas.ipcc.ch/regional-information) to ensure that the GCMs used to develop the multi-model ensemble encompassed a range of ocean temperature futures projected for eastern Australia under climate change. CMIP6 SST projections under the SSP5-8.5 scenarios were extracted from the CMIP6 data portal for analysis (see https://esgf-node.llnl.gov/search/cmip6).

The coarse spatial resolution of GCM data (0.25–1°) has limited utility for projecting changes in suitable thermal habitat for each study species at popular dive locations along the eastern Australian coastline. Therefore, we applied the delta ‘change-factor’ method to downscale future SST data throughout the study region to a 0.05° spatial resolution that is consistent with satellite SST observations. Delta downscaling was undertaken because it has proven utility for providing high-resolution mean climate conditions over decadal time periods for numerous climate-impact studies (e.g. Morley et al. 2018; Navarro-Racines et al. 2020; von Hammerstein et al. 2022).

This downscaling process involved the following: (1) re-mapping the curvilinear source GCM data to a global rectilinear grid by using the second-order conservative algorithm (remapcon2) in Climate Data Operators (ver. 1.9.6, see https://doi.org/10.5281/zenodo.2558193; Schulzweida 2019); (2) calculating the difference (i.e. delta value) between monthly aggregated data for successive 20-year periods centred on all years ranging between 2015 and 2059 and a modelled historical baseline period encompassing 1993–2012 for each GCM; (3) disaggregating delta value matrices from their native model resolution to the finer resolution of satellite-observed SST data (i.e. 0.05°) by using bilinear interpolation; and (4) adding delta value matrices centred on years ranging between 2015 and 2059 to an observed monthly SST climatology encompassing the period 1993–2012. Satellite SST data used consistently throughout this study are daily global SST from the CMEMS (product #010_011). A multi-model ensemble was then created using median SST values from downscaled projections across each of the seven GCMs. Median values were utilised to create the multi-model ensemble because this statistic is robust to outliers and does not assume that the data are normally distributed. This method produced a monthly time-series (2015–59) of future SST data downscaled to a 0.05° resolution required to facilitate projections of thermal-habitat suitability for manta rays and zebra sharks off eastern Australia.

To assess the accuracy of future habitat projections produced using the ensemble for CMIP6 SST data, monthly thermal-habitat suitability matrices nearshore of the continental shelf break were generated using satellite-derived and modelled SST for the 6-year period encompassing 2015–20. Collinearity among monthly spatial matrices of species thermal-habitat suitability created using observed and modelled data for the corresponding historical period (i.e. 2015–20) was assessed using Pearson’s correlation coefficient to infer the accuracy of future projections of thermal-habitat suitability (Spillman and Alves 2009; Eveson et al. 2015).

Spatial matrices of thermal-habitat suitability, for both study species, created nearshore of the continent shelf break at monthly timesteps between January 2015 and December 2020 by using observed and modelled SST data sources were positively correlated (Supplementary Fig. S1). Pearson correlation coefficients evaluating the collinearity of thermal-habitat suitability matrices produced using observed and modelled SST data ranged between 0.77 and 0.97 for manta rays and 0.47–0.95 for zebra sharks. This supports the use of thermal-habitat suitability projections created using downscaled SST data from an ensemble of CMIP6 GCMs for assessing future changes in manta ray and zebra shark temporal persistence off eastern Australia.

Temporal persistence at dive locations

Monthly spatial projections of manta ray and zebra shark thermal-habitat suitability, from January 2015 to December 2059, were used to assess changes in annual temporal persistence (i.e. number of months per year) at popular dive locations in south-eastern Australia. The dive locations assessed were Byron Bay (28.64°S), Coffs Harbour (30.30°S), South West Rocks (30.92°S), Forster (32.19°S), Port Stephens (32.69°S), Sydney (33.89°S) and Jervis Bay (35.19°S). These locations were selected because they are popular scuba-diving locations in NSW with currently operating recreational scuba-diving businesses (Harriott et al. 1997; Willis and Boshoff 2020). For each monthly projection, thermal-habitat suitability at each dive location had to exceed a species threshold that distinguishes suitable from unsuitable thermal habitat. Given that habitat models predicted thermal suitability on a scale ranging between 0 and 1, it was necessary to identify a threshold value for partitioning predictions into ‘suitable’ and ‘unsuitable’ habitat categories. A mid-point threshold of 0.5 is commonly used for this purpose in habitat-modelling studies (Bailey et al. 2002; Stockwell and Peterson 2002); however, arbitrarily adopting this value has been criticised for lacking ecological meaning (Liu et al. 2005). Instead we created daily predictions of thermal-habitat suitability for each study species at locations (n = 50) of their known occurrence randomly sampled from the historical occurrence dataset. The median of these values was then selected as the threshold for suitable thermal conditions for each species. Median values were selected (as opposed to minimum values or low quartiles) because these provide an estimate of the thermal conditions capable of supporting the occurrence of high abundances of the study species, rather than a lower level of thermal-habitat suitability that may support only the occurrence of smaller numbers of individuals in space and time. This was considered conservative, given that our analyses aim to relate changes in the temporal persistence of suitable thermal conditions for the study species to opportunities that divers have to encounter these species, and sufficiently high abundances are required to facilitate encounters. Furthermore, applying this data-driven approach for identifying threshold values maximises the agreement between modelled and observed species distributions, ultimately optimising the accuracy of temporal-persistence estimates (Champion et al. 2019).

Changes in the temporal persistence of suitable thermal conditions for manta rays and zebra sharks at dive locations over the future projections period was assessed using linear models. Slope coefficients and confidence intervals from model fits were used to estimate rates of change in the number of months per year that suitable thermal conditions were present at dive locations. P-values were assessed to determine whether slopes for each species at each dive location were significantly different from zero. Residuals from these models were linearly related to theoretical quantiles and no patterns in residual variation were evident, indicating that our analyses satisfied statistical assumptions. Data analysis was undertaken using the R software (ver. 4.2.3, R Foundation for Statistical Computing, Vienna, Austria, see https://www.r-project.org/), with GAMMs fitted using the gamm4 package (ver. 0.2-6, see https://cran.r-project.org/package=gamm4), k-fold cross-validation undertaken using the dismo package (ver 1.3-9, see https://cran.r-project.org/package=dismo), spatial data aggregated using the raster package (ver. 3.6-20, see https://cran.r-project.org/package=raster) and figures generated using the ggplot2 package (ver. 3.4.1, see https://CRAN.R-project.org/package=ggplot2; Wickham 2016).

Habitat-suitability models

Multiple dynamic oceanographic variables were found to be significant linear and non-linear predictors of the occurrence of manta rays and zebra sharks off eastern Australia (Fig. 2, Table 3). SST had a unimodal effect on environmental habitat suitability for both study species, with peaks in optimal thermal habitat occurring between 23 and 24.5°C for manta rays and 24.5 and 26°C for zebra sharks (Fig. 2). The occurrence of both study species was associated with a negative linear response to EKE (i.e. current velocity), indicating a preference for environments characterised by low to moderate current velocity. Zebra shark occurrence was also associated with a negative linear response to SLA, whereas there was no significant association between manta ray occurrence and SLA. Both species also responded to the concentration of chlorophyll in seawater, with increases in chlorophyll concentration being generally associated with non-linear improvements in environmental habitat suitability (Fig. 2).

Fig. 2.

Partial effects of covariates on the fitted values of the best manta ray (Mobula alfredi) and zebra shark (Stegostoma tigrinum) environmental-suitability models. Grey shading denotes 95% confidence intervals around the model fit (black line). SST, sea-surface temperature; EKE, eddy kinetic energy; SLA, sea-level anomaly; CHL, chlorophyll-a concentration.

Importantly, SST was found to be a dominant predictor of the occurrence of manta rays and zebra sharks off eastern Australia. Six-fold model cross-validation showed that only minor reductions in model predictive skill, as assessed using AUC (Δ < 0.07) and TSS (Δ < 0.06) measures, occurred when the best models were reduced to temperature-only models (Table 4). Mean AUC and TSS values associated with temperature-only models for both species consistently exceeded 0.75 and 0.45 respectively (Table 4). These values denote models associated with fair to good performance (Swets 1988; Pearce and Ferrier 2000), indicating that our future projections of species thermal-habitat suitability are likely to provide accurate approximations of species spatial distributions under ocean warming.

Spatial projections of thermal-habitat suitability

Spatial analyses highlighted that suitable thermal habitat for manta rays and zebra sharks undergo a poleward shift to higher latitudes along Australia’s eastern coast during the Austral summer and autumn, and an equatorward retreat during the cooler winter and spring seasons (Fig. 3, 4). Seasonal shifts in the distribution of thermal habitat for manta rays and zebra sharks along the eastern coast of Australia are projected to be exacerbated under future ocean warming, with suitable thermal habitat for both species projected to occur further south in all seasons under ocean warming (Fig. 3, 4).

Fig. 3.

Thermal-habitat suitability for manta rays (Mobula alfredi) off eastern Australia predicted for (a) the 20-year period centred on 2010, and projected for (b) the 20-year period centred on 2050. (c) The change in thermal-habitat suitability for manta rays between 2010- and 2050-centred periods. Monthly spatial predictions were averaged for each period and seasonally aggregated (summer, December–February; autumn, March–May; winter, June–August; spring, September–November). Predictions have been made nearshore of the continental shelf-break (i.e. 200-m isobath).

Fig. 4.

Thermal-habitat suitability for zebra sharks (Stegostoma tigrinum) off eastern Australia predicted for (a) the 20-year period centred on 2010, and projected for (b) the 20-year period centred on 2050. (c) The change in thermal habitat suitability for zebra sharks between 2010- and 2050-centred periods. Monthly spatial predictions were averaged for each period and seasonally aggregated (summer, December–February; autumn, March–May; winter, June–August; spring, September–November). Predictions have been made nearshore of the continental shelf-break (i.e. 200-m isobath).

Future temporal persistence at dive locations

Significant increases in the temporal persistence of suitable thermal habitat for manta rays and zebra sharks are projected at all dive locations distributed between Forster (32.19°S) and Jervis Bay (35.19°S) under future ocean warming (Fig. 5, Table 5). Furthermore, a significant increase in the temporal persistence of suitable thermal habitat for zebra sharks at Coffs Harbour is also projected in the future (Fig. 5, Table 5). Contrastingly, no significant changes in suitable thermal habitat for manta rays are projected for sites from South West Rocks to Byron Bay and for zebra sharks at South West Rocks and at Byron Bay (Fig. 5, Table 5).

Fig. 5.

Temporal trends in the projected annual persistence of suitable thermal habitat for manta rays (left-side panel) and zebra sharks (right-side panel) at key dive locations from Australia’s eastern coast. Data presented are based on monthly spatial projections of thermal-habitat suitability by using an ensemble of downscaled sea-surface temperature data forced under the SSP5-8.5 scenario between 2015 and 2059.

The number of months per year that suitable thermal conditions for manta rays are projected to occur, are projected to increase between 2015 and 2059, with increases ranging between ~2.5 (i.e. at Jervis Bay) and 4 (i.e. at Sydney) months per year. For zebra sharks, increases of between ~0.5 and 3 months per year are projected across dive locations in the future-projection period (Table 5). Importantly, the absolute number of months per year that suitable thermal conditions are projected to occur in the future follows a latitudinal trend for both study species, with larger increases in the number of months per year projected at lower-latitude dive locations, albeit with this trend being more consistent for manta rays (Table 5, Fig. 5).

Effect of environmental variables on current species-migration patterns

Sea-surface temperatures were found to be the dominant factor associated with the presence of manta rays and zebra sharks off eastern Australia, with suitability for the occurrence of manta rays peaking when sea-surface temperatures range between 23 and 24.5°C for manta rays and between 24.5 and 26°C for zebra sharks. Spatial analyses highlighted that suitable thermal conditions for these species undergo a poleward shift towards higher latitudes during the austral summer and autumn, and an equatorward retreat during the cooler winter and spring seasons. Our results are consistent with previous research on movements of manta rays, conducted at Lady Elliott Island in Queensland, which found that some manta rays migrate south, to North Stradbroke Island and Byron Bay, from mid-spring to mid-autumn (Couturier et al. 2011). Model predictions are also consistent with seasonal occurrence data for zebra sharks in south-eastern Queensland, which identified that zebra sharks start to move into this area in November, once ocean temperatures reached 22°C and then migrate further south as temperature increased above this threshold (Dudgeon et al. 2013).

In addition to moving in response to shifting ocean temperatures, occurrences of manta rays and zebra sharks were also found to vary with current velocity and chlorophyll concentration, with both species exhibiting a preference for regions with low to moderate current velocities and increased chlorophyll concentrations. Increased presence of manta rays in areas with high chlorophyll concentrations is understandable, given that high chlorophyll concentrations generally indicate high concentrations of plankton (Jakobsen and Markager 2016) and manta rays are planktonic filter feeders that can migrate large distances to feed in plankton-rich waters (Graham et al. 2012; Jaine et al. 2014). For zebra sharks, increased presence in regions with high chlorophyll concentrations is harder to explain, although potentially higher chlorophyll concentrations and associated high plankton concentrations may increase abundance of their food (molluscs), with molluscs forming the bulk of the diet for zebra sharks (Cortés 1999). The effects of current strength on manta ray and zebra shark distributions are likely to be more complex, although it can be argued that light–moderate currents provide preferrable conditions, allowing easy movement and feeding, whereas stronger currents are known to inhibit fish movement (Caldwell and Gergel 2013), with strong currents also being associated with lower manta ray sightings in Mozambique (Rohner et al. 2013).

Implications for fisheries management

In addition to affecting dive tourism in NSW, changes to the migration patterns for manta rays and zebra sharks have implications for fisheries management. Zebra sharks are known to be caught as bycatch in Queensland (Chin et al. 2012) and elsewhere (Pottie et al. 2021). Manta rays have been recorded as bycatch under the NSW Shark Meshing (Bather Protection) Program (SMP; Broadhurst and Cullis 2020), with 0.30 entanglements per year from 1950 to 1993 (Krogh and Reid 1996) and 0.53 per year from 2009 to 2024 (M. Green, pers. comm.). Given the projected increased presence and southward range extensions of these species in NSW, there is potential for increasing exposure to SMP and fishing gears. Consequently, fisheries management plans should continue to monitor interactions between species and SMP and other fishing gear into the future and examine whether changes to management policies are needed to address these and any other identified threats.