Future Distribution of Suitable Habitat for Pelagic Sharks in Australia Under Climate Change Models

Global oceans are absorbing over 90% of the heat trapped in our atmosphere due to accumulated anthropogenic greenhouse gases, resulting in increasing ocean temperatures. Such changes may influence marine ectotherms, such as sharks, as their body temperature concurrently increases toward their upper thermal limits. Sharks are high trophic level predators that play a key role in the regulation of ecosystem structure and health. Because many sharks are already threatened, it is especially important to understand the impact of climate change on these species. We used shark occurrence records collected by commercial fisheries within the Australian continental Exclusive Economic Zone (EEZ) to predict changes in future (2050–2099) relative to current (1956–2005) habitat suitability for pelagic sharks based on an ensemble of climate models and emission scenarios. Our predictive models indicate that future sea temperatures are likely to shift the location of suitable shark habitat within the Australian EEZ. On average, suitable habitat is predicted to decrease within the EEZ for requiem and increase for mackerel sharks, however, the direction and severity of change was highly influenced by the choice of climate model. Our results indicate the need to consider climate change scenarios as part of future shark management and suggest that more broad-scale studies are needed for these pelagic species.


INTRODUCTION
Climate change is predicted to have unprecedented effects on the marine environment, with changes in ocean temperature increasing extinction risk for many species (Dulvy et al., 2003;Barnosky et al., 2011;Bruno et al., 2018;Pinsky et al., 2019) and altering the global distribution of marine life (Tittensor et al., 2010;Garciá Molinos et al., 2016). Changes in species distribution (Perry et al., 2005;Poloczanska et al., 2013) and community structure (Doney et al., 2012) are already being observed in marine ecosystems due to temperature shifts associated with rising emissions and accumulation of atmospheric carbon dioxide (Hoegh-Guldberg and Bruno, 2010;Doney et al., 2012;Gattuso et al., 2015). Recent modeling of biodiversity under different future climate change scenarios, across a wide range of marine and terrestrial ecosystems, predicts abrupt and irreversible ecosystem disruption during the late 21st century (Trisos et al., 2020). With predicted increases of up to ∼5 • C in worldwide sea-surface temperature (SST) by the end of the 21st century (IPCC, 2015), there is a critical need to investigate how marine species will be affected, especially ectotherms which are dependent on external sources for body heat. As ectotherms, sharks may be influenced by climate change (Bernal et al., 2012;Rosa et al., 2014Rosa et al., , 2017Syndeman et al., 2015;Pinsky et al., 2019), with higher temperatures increasing their metabolism and oxygen demand (Pistevos et al., 2015;Lawson et al., 2019). The exception to this may be Lamnid mackerel sharks, which have some endothermic capability (Watanabe et al., 2015).
Many shark species are already globally threatened due to fisheries overexploitation (Queiroz et al., 2019) coupled with their low fecundity, late age at maturity, and slow growth (Cortés, 2000;Garcia et al., 2008;Yokoi et al., 2017). In fact, 16.6% of shark species are estimated to be threatened with extinction, and another 37.9% of shark species are categorized as "Data Deficient" by the International Union for Conservation of Nature (IUCN, 2020). Nevertheless, sharks are known to have direct economic value in fisheries (Dulvy et al., 2017) and ecotourism (Cisneros-Montemayor et al., 2013;Huveneers et al., 2017). They also play a key role in ecosystem functioning and stability, connecting distant ecosystems via their long-distance migrations (Rogers et al., 2015), and altering prey behavior, distribution and energy use (Heupel et al., 2015;Roff et al., 2016;Dulvy et al., 2017). Climate change may exacerbate existing threats for sharks, for example, suitable pelagic shark habitat in the north Pacific Ocean is projected to decline by the year 2100 (Hazen et al., 2013).
Future projections based on existing observations and modeling techniques can be used to investigate the effects of climate change on pelagic sharks (Barange et al., 2016). Using Earth System Models from the Coupled Model Intercomparison Project Phase 5 (CMIP5; hereafter called "climate models"), complex relationships between ecosystem health, human activities and global climate can be included to evaluate alternative future scenarios with varying severity of emissions (Moss et al., 2010;Freer et al., 2017). There are four emission scenarios commonly referred to as Representative Concentration Pathways (RCP 2.6, RCP 4.5, RCP 6, and RCP 8.5) (IPCC, 2013). These RCP scenarios are used to predict radiative forcing values, a measure of absorbed and retained energy in the lower atmosphere, for the year 2100 (Moss et al., 2010;Vuuren et al., 2011). RCP 4.5, also referred to as "stabilization scenario, " is an optimistic scenario assuming a decline in overall energy usage from fossil fuel sources that limits emissions and radiative forcing . Conversely, RCP 8.5, also referred to as "business-as-usual, " is the most pessimistic scenario assuming minimal stabilization of greenhouse gas emissions alongside a large human population with high energy demands .
The Australian Exclusive Economic Zone (EEZ) is already being impacted by climate change with waters off south-east Australia warming at almost four times the global average (Oliver et al., 2017) and range extensions already documented in several fish species (Last et al., 2011). Australia has one of the world's most diverse communities of sharks, with 182 recognized species (Simpfendorfer et al., 2019), and SST has been shown consistently to be a strong predictor of pelagic shark occurrence in Australian waters (Rogers et al., 2009(Rogers et al., , 2015Stevens et al., 2010;Heard et al., 2017;Birkmanis et al., 2020). It is therefore important to investigate the likely impact of temperature changes on pelagic shark distribution and the location of suitable habitat on a continental scale if these species are to be appropriately managed into the future -especially if such changes may require a reassessment of interactions with fisheries in the future. Sharks comprise approximately 27% of the total catch (by number) of Australian pelagic longline fisheries (Gilman et al., 2008), with Australian stocks of the IUCN classified "Critically Endangered" oceanic whitetip (Carcharhinus longimanus), "Endangered" shortfin mako (Isurus oxyrinchus), and "Endangered" longfin mako (Isurus paucus) sharks listed respectively as "overfished, " "depleting, " and "undefined" due to a lack of data (Simpfendorfer et al., 2019;IUCN, 2020).
This study follows on from Birkmanis et al. (2020) in which occurrence records of pelagic sharks belonging to the Carcharhinidae and Lamnidae families (hereafter "requiem" and "mackerel, " respectively) were obtained from commercial fisheries and used to develop generalized linear models with which to predict suitable habitat for these species within the Australian continental EEZ. After accounting for fishing effort bias, these models showed that SST was an important predictor of shark distributions, with the highest ranked model also including turbidity. Here, we extend our modeling to assess the impact of climate change on pelagic shark habitat across the entire continental Australian EEZ.

Shark Occurrence
Catch records of 3,973 individual sharks from two families; requiem (silky Carcharhinus falciformis, oceanic whitetip Carcharhinus longimanus, dusky Carcharhinus obscurus and blue Prionace glauca) and mackerel (shortfin mako Isurus oxyrinchus, longfin mako Isurus paucus and porbeagle Lamna nasus) were obtained through the Global Biodiversity Information Facility online database (GBIF.org, 2017), as per details included in Birkmanis et al. (2020). These oceanic sharks were caught predominantly using commercial longlines in Commonwealth managed fisheries (more detailed data unavailable), with catch locations depicted in Supplementary Figure S1.

Predictors for Modeling Baseline and Future Climate Environmental Data
A climatological baseline was used as a reference point for projected future climate changes. According to Birkmanis et al. (2020), SST and turbidity were the most suitable predictors of requiem and mackerel shark occurrence within the Australian EEZ. We therefore focused on these two predictors to develop a climatological baseline to use as a reference for projected future climate changes. To calculate the SST baseline data, we downloaded monthly SST values for the years 1956-2005, covering the time period of our observed shark occurrence data, from the Integrated Marine Observing System (IMOS, 2016). We then averaged the SST values for each 0.1 • gridcell in the study area using ArcGIS 10.5 from Environmental Systems Research Institute (ESRI, 2017). We incorporated the observed turbidity values (measured as mean diffuse attenuation coefficient at wavelength 490 nm, downloaded using the Marine Geospatial Ecology Tool; Roberts et al., 2010) from 2000 to 2002 into our models with the assumption that turbidity will remain unchanged in the future.
Future SST data were taken from 24 CMIP5 climate models, using only one realization per climate model, under two emission scenarios, RCP 4.5 and RCP 8.5, amounting to 48 total simulations ( Table 1). We downloaded the SST field and the anomaly statistic for each climate model (  Scott et al., 2016 for details). We used the portal to calculate the difference in the mean SST between the future climate (2050-2099) and the model baseline reference period , hereafter called "anomaly" data. We then added these anomaly data to our baseline data across the extent of the Australian EEZ using ArcGIS and included this as the SST predictor for the future values.

Modeling Habitat Suitability for Baseline and Future Climate Data
We developed binomial generalized linear models with a logit link function for each of the two pelagic shark families, following  Birkmanis et al. (2020). In brief, the probability of shark occurrence (calculated as the number of sharks caught divided by the number of fishing boats occurring in the same gridcell) was used as the response variable, with turbidity and SST values for either the climatological baseline or the future used as predictors. We included effort, defined as the number of boats recorded in each grid-cell from the same time period as the occurrence data (2000)(2001)(2002), as a model weight to account for differing amounts of catch per unit effort (CPUE) within the entire EEZ. As in, we weighted our models by fishing effort to estimate the probability of finding a shark in each grid-cell within the Australian EEZ which minimized the effect of fisheries effort on the data (Birkmanis et al., 2020). To stabilize parameter estimation, we standardized both predictors to z-scores using the scale function in R statistical software (R Core Team, 2017) before inclusion in our models (James et al., 2015). We also included a quadratic term for SST using the poly function from the stats package (R Core Team, 2017) in R statistical software to account for likely preferential SST ranges. We then quantified the goodness-of-fit for each model using the percentage of deviance explained, and used the predict function from the stats package in R statistical software to predict shark habitat suitability for the baseline data and also for the end of the century using the future climate data. To calculate the amount of change in suitable habitat under each climate model and emission scenario, we subtracted the number of grid-cells with resulting suitable habitat ≥0.5 in the future climate scenarios from those obtained in the baseline scenario. Differences between baseline and future scenarios show the change in suitable habitat area predicted for each family under possible future conditions.

RESULTS
Anomalies in SST in the Australian EEZ varied according to the climate model and emission scenario used (Figure 1 and Table 1).
The mean SST anomaly for all climate models was 2.27 • C (SD: 0-1.2) for RCP 4.5 and 3.78 • C (SD: 0-1.21) for RCP 8.5. Our results show that the predicted mean SST anomaly ranged from minima of 0.93 • C (for climate model GISS-E2-R, RCP 4.5) and 0.69 • C (for climate model HADGEM2-AO, RCP 8.5), to mean maxima of 1.83 • C (for climate model IPSL-CM5A-MR for both emission scenarios) at the end of the century (Figure 1 and Table 1).
The climate model MPI-ESM_LR resulted in the maximum SST anomaly projected by all climate models (3.76 • C for RCP 4.5 and 5.71 • C for RCP 8.5, respectively). Despite model-to-model variation in the magnitude of anomalies, all climate models predicted south-eastern Australia would experience the greatest SST increases by the end of the century (Figure 1). Both the baseline and future habitat suitability models explained slightly higher deviance for requiem than mackerel sharks but all values were around 30% (Supplementary Table S1). The baseline models explained 31.13 and 27.33% for requiem and mackerel sharks, respectively. The future models explained 29.91 and 31.76% for RCP 4.5 and RCP 8.5, respectively for requiem, and 26.47 and 26.22% for RCP 4.5 and RCP 8.5, respectively for mackerel sharks. The resulting predicted habitat suitability maps are presented as the mean across all climate models for requiem (Figure 2) and mackerel sharks (Figure 3), with the predicted change per climate model presented in  requiem and mackerel sharks (0.65 and 0.63, respectively). For both requiem and mackerel sharks, the maximum habitat suitability (∼0.8) was predicted by climate model NORESM1-ME under both emission scenarios (Figure 4). Regions where habitat was predicted to be suitable (i.e., ≥0.5) at the end of the century varied by family, with southern Australia suitable for mackerel sharks, and north-eastern Australia for requiem sharks (Figures 2, 3).
Based on 48 climate simulations, our results suggest a shift in suitable habitat for both requiem and mackerel sharks within the Australian EEZ in the last half of the twenty-first century . The severity and direction of this shift varied, with suitable habitat for requiem sharks predicted to decrease under most climate models, while habitat suitability for mackerel sharks varied to a greater degree depending on the climate model and emission scenario. On average, predicted suitable habitat for requiem sharks under RCP 4.5 extended south on the north-eastern (∼600 km) and south-western coast (∼200 km), but decreased in the north-west (∼400 km). For RCP 8.5, suitable habitat was projected to extend south on the northeastern coast (∼650 km) and decrease across the north-west (∼500 km) with similar increases on the south-western coast (Figure 2). For mackerel sharks, the average of all climate models predicted an increase in suitable habitat across on the southern coast (∼900 km) and off the southern extent of the EEZ south of Tasmania (∼400 km) for RCP 4.5, with increases also projected to occur under RCP 8.5 (∼700 km across and ∼200 km south along the southern coast and ∼150 km south off the southern extent of the EEZ south of Tasmania) (Figures 3, 4).

DISCUSSION
Significant shifts in the distributions of marine organisms are being observed in the global ocean due to anthropogenic climate change (Poloczanska et al., 2013). Our results highlight that shifts in the location of suitable habitat for requiem and mackerel sharks by the end of the century are to be expected, with a decrease in predicted suitable habitat for requiem sharks off the south-western coast under both emission scenarios. This agrees with predicted habitat shifts for silky, blue (both requiem family; Cheung et al., 2015;Lezama-Ochoa et al., 2016), and mako sharks (mackerel family; Hazen et al., 2013) in other areas. The waters of south-western and south-eastern Australia are warming at an increased rate, almost three and four times higher than the global average, respectively (Hartmann et al., 2013;Robinson et al., 2015a) as indicated in Figure 1. Our models predict that this area will become unsuitable for both requiem and mackerel sharks, likely due to the water temperatures at the end of the century exceeding the thermal tolerance of these pelagic sharks. In our analysis, and those of Robinson et al. (2015b) and Hobday (2010), southward shifts in suitable habitat for blue and mako sharks on the eastern coast of Australia are predicted. This is in line with ocean climate zones (areas with distinct climate, based on annual SST values) shifting southwards by 200 km along the north-eastern coast and approximately 100 km along the northwestern coast in tropical Australian waters (Lough, 2008). In the north Pacific Ocean suitable habitat loss was predicted for both blue and mako sharks by the end of the century (Hazen et al., 2013). Such differences in predictions may be due to currents and northern latitude prey species being able to migrate poleward along the coastline (Perry et al., 2005). Due to the east-west orientation of the temperate Australian coastline and limited continental shelf area to the south of the continent (Urban, 2015), there are few opportunities for continental shelf marine organisms, including fish that are shark prey species, to move to higher latitudes and avoid increased water temperatures. Even with suitable habitat available for pelagic sharks within Australian waters these predators will follow prey species, such as tuna (Hobday and Poloczanska, 2010), which are expected to decline in the tropics and shift poleward in response to a warming ocean (Erauskin-Extramiana et al., 2019).
Although relatively little is known about how elevated temperatures will affect sharks (Pistevos et al., 2015), pelagic sharks are vulnerable to climate change impacts (Jones and Cheung, 2018) and life history strategies may play a part in determining ultimate patterns of species distribution. For relatively sedentary, benthic shark species, exposure to projected end-of-century temperatures has been shown to result in both positive and negative impacts. Port Jackson sharks (Heterodontus portusjacksoni) exposed to elevated temperatures exhibited an increase in mortality, altered behavior, increased learning performance and feeding, but reduced growth and embryonic development time (Pistevos et al., 2015;Vila et al., 2018Vila et al., , 2019. Conversely, brownbanded bamboo sharks (Chiloscyllium punctatum) showed decreased survival alongside significantly increased embryonic growth and ventilation rates (Rosa et al., 2014), while juvenile epaulet sharks (Hemiscyllium ocellatum) showed significantly decreased growth rates and 100% mortality . It is likely that the physiological impacts of increasing ocean temperature will be greater for more active pelagic sharks than for benthic species (Rosa et al., 2014), given their reliance on ram ventilation and continuous movement (Lawson et al., 2019). Sharks already at their provisioning limit may be faced with starvation if temperature-driven FIGURE 4 | Change in predicted habitat suitability for requiem and mackerel sharks in the Australian EEZ between the baseline time period  and at the end of the twenty-first century (2050-2099) under two emission scenarios (RCP 4.5 and 8.5). increases in metabolic rates are not met with higher food intake (Pistevos et al., 2015), and this risk will be heightened should environmental perturbations concurrently influence prey availability and abundance. However, the thermal tolerance of requiem and mackerel sharks (Francis and Stevens, 2002;Last and Stevens, 2009;Corrigan et al., 2018;Hueter et al., 2018;Young and Carlson, 2020) may enable them to cope with changing temperatures.
Even though we predicted an overall increase in the amount of suitable habitat for mackerel sharks at the end of the 21st century, temperature acclimatization comes with an energetic cost that impacts other functions such as reproduction, growth, foraging and swimming. Changes in the marine environment may result in novel ecosystems requiring predators to alter foraging behaviors and adapt to new prey species (Nagelkerken and Munday, 2016;Rivest et al., 2019). Under such stresses, individuals become less competitive with decreases in reproduction and population density (Beaugrand and Kirby, 2018) and may exploit habitat heterogeneity by undertaking vertical migrations to suitable temperatures to maximize biological efficiency and minimize physiological adjustment costs (Chin et al., 2010;Beaugrand and Kirby, 2018). The endothermic ability to swim faster and farther (Watanabe et al., 2015) may allow mackerel sharks to migrate longer distances and forage over wider areas with greater access to prey and seasonal resources, although at higher energetic costs than ectothermic species. However, the ability of pelagic sharks to move and follow the shifting suitable habitat outside their current ranges, may potentially alter their interactions with fisheries. It is worth noting that latitudinal species shifts in response to warming can be misleading with some pelagic species migrating vertically not latitudinally (Perry et al., 2005;Beaugrand and Kirby, 2018) and this may be the case with some pelagic shark species. In Australian waters, pelagic sharks have been recorded regulating their depth to occupy regions of favorable temperatures, although this behavior could also be related to prey movements (Rogers et al., 2009;Stevens et al., 2010;Heard et al., 2017) as well as habitat suitability.
Our study predicts changes in habitat suitability for pelagic sharks in the Australian EEZ, but predictions at the end of the century are highly dependent on the climate model and emission scenario chosen to represent future conditions. This has been the case for similar studies on other species, for example, freshwater fish assemblages (Buisson et al., 2010) and mesopelagic lanternfish in the Southern Ocean (Freer et al., 2019) highlighting the benefit of using an ensemble approach to capture high climate uncertainty. Moreover, the SST anomalies across the Australian EEZ also vary according to the climate model and emission scenario used in the analysis. Our analysis was done at the family level due to the sample size available. Analysis at family level, whilst valuable for relatively homogenous species groups, inevitably results in loss of information at lower taxonomic levels. Further research is needed in more localized areas, including telemetry studies on single species, to add greater certainty to species distribution model predictions. There is no consensus about how turbidity may vary under a changing climate, and in our models we assumed that turbidity levels would remain stable at the end of the century. However, a predicted increase in extreme rainfall influenced by changes in atmospheric circulation may increase coastal turbidity due to terrestrial-derived nutrients and pollutant input (Harley et al., 2006). Additionally, turbidity is correlated with chlorophylla in pelagic systems, and warmer water temperatures drive phytoplankton blooms, with elevated temperatures increasing both cyanobacterial and algal chlorophyll-a concentrations (Lürling et al., 2018;Trombetta et al., 2019). As aquatic nutrients have a greater impact on chlorophyll-a concentrations than temperature, and salinity and wind are also correlated with plankton blooms (Lürling et al., 2018;Trombetta et al., 2019), the impact this may have on pelagic systems in Australian waters is still unclear. Despite the uncertainties associated with predicting future conditions, studies such as ours using remotely sensed environmental information and occurrence data from fisheries over a large spatial scale, are important to understand how pelagic species with broad geographic ranges might fare in the future. Such studies are a first component of broader research in which the distribution of multiple species are predicted in a likely altered future marine environment.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
CB, AS, and JP conceived the study. CB, JP, AS, and LS designed the methodology with assistance from JF. CB collated and analyzed the data and led the writing of the manuscript with significant contributions from AS. All authors contributed critically to the drafts and gave final approval for publication.

ACKNOWLEDGMENTS
We thank R. Summerson for assistance with accessing the fisheries data and acknowledge Australian Bureau of Agricultural and Resource Economics and Sciences (ABARES) as the source of the fisheries data, originally supplied by Australian Fisheries Management Authority (AFMA) and state fisheries management agencies.