Seascape topography slows predicted range shifts in fish under climate change

Shifts in marine species distributions are occurring rapidly in response to climate change, yet predicted rates of change are partly dependent on the data used to estimate species distributions. While pelagic fishes are known to respond to dynamic oceanographic variables, static topographic features can also regulate their distributions, yet the effect of topography on rates of range shifts remains poorly understood. For a suite of five pelagic fishes from an ocean warming hotspot, we developed generalized additive mixed models that did and did not incorporate preferences for seascape topographic heterogeneity when estimating species distributions. Monthly spatial predictions over two decades revealed that poleward extensions in species' core environmental habitats were significantly reduced by 5.3–60.4% (7.1–131.4 km per decade) when preferences for seascape topography were incorporated. These findings suggest that ignoring species associations with static environmental features can substantially affect predicted rates of climate‐driven redistributions in fishes.

Shifts in species geographic distributions have emerged as one of the most globally evident biological responses to climate change, with far reaching implications for biodiversity conservation and human well-being (Pecl et al. 2017). Rates of climate-driven redistributions suggest that marine species are shifting toward high latitudes between 6 and 10 times more rapidly than terrestrial taxa (Lenoir et al. 2020). However, considerable variability exists among rates of range shifts in marine systems (Poloczanska et al. 2016), with regional rates of environmental warming accounting for a greater proportion of variation in species responses than biological traits (Pinsky et al. 2013).
Species distribution models (SDMs) are value tools for quantifying range shifts in marine environments, particularly for highly mobile species (Robinson et al. 2017). Range shift analyses utilizing SDMs commonly focus on quantifying distributional shifts in dynamic oceanographic habitats (Champion et al. 2018;Hobday 2010) that are determined by species' preferences for variables such as sea surface temperature, mixed layer thickness and sea surface height. However, range shift analyses based on dynamic oceanographic habitats do not consider species' associations with static environmental habitats, such as seascape topography, which are also known to structure distributions and abundances of marine species (Hobday and Campbell 2009;Morley et al. 2018;Mejía-Mercado et al. 2019). Given that the selection of covariates used to capture species environmental habitat associations can significantly affect model predictions (McHenry et al. 2019;Brodie et al. 2020), range shift rates estimated using models that exclude static environmental variables may not be representative of species realized distributions. Improving the accuracy of range shift predictions for marine species requires quantitative comparisons that account for both dynamic oceanographic processes and static environmental habitat associations.
Pelagic fishes are undergoing some of the most rapid climate-driven redistributions reported for any group of marine species globally (Hazen et al. 2013;Poloczanska et al. 2013). For example, rates of recent poleward redistributions of pelagic fishes off eastern Australia (e.g., Seriola lalandi: 94.4 km per decade [Champion et al. 2018] and Istiompax indica: 88.2 km per decade [Hill et al. 2015]) markedly exceed the average rate of range shifts in nearshore fishes from the same region (38 km per decade; Sunday et al. 2015). Given that dynamic oceanographic variables are good predictors of the distribution and body condition of pelagic fishes (Zainuddin et al. 2006;Champion et al. 2020), range shift analyses for pelagic fishes have focused on spatiotemporal shifts in species' dynamic oceanographic habitats (Hill et al. 2015). However, pelagic fishes also associate with topographic features such as reefs and seamounts despite being highly mobile and wide-ranging (Hobday and Campbell 2009;Rees et al. 2018). For example, catch rates of pelagic fishes off southern Australia were up to 26 times greater near topographic features than away from them (Hobday and Campbell 2009), demonstrating a strong influence of topographic heterogeneity on pelagic fish distributions. Explanations for topographic habitat associations in pelagic fishes include refuge from predation, formation of spawning and feeding aggregations and cleaning stations, and also the provision of waypoints on migration routes (Biesinger et al. 2011;Fréon and Dagorn 2000). Given that continental shelf seascapes commonly contain topographic features that attract pelagic fishes and potentially alter rates of redistribution, omitting topographic habitat associations from range shift analyses may substantially affect predictions of climate-driven range shifts.
We propose that topographic habitat associations mediate rates of predicted climate-driven range shifts in pelagic fishes. Using five pelagic fishes from an ocean warming hotspot, we quantified species' topographic habitat preferences and developed SDMs that both omitted and included these environmental responses to measure the relative effect of topography on rates of climate-driven redistributions. Our findings highlight the importance of incorporating static (e.g., topographic) variables into analyses predicting the climate-driven redistribution of pelagic fishes, and have implications for assessing the vulnerability of pelagic fishes to changing ocean conditions.

Case-study system and species
The spatial extent of this study encompassed the marine environment adjacent to eastern Australia (145-160 E; 15-45 S; Fig. 1). Climate-driven strengthening of the East Australian Current is causing the penetration of anomalously warm seawater to higher latitudes off eastern Australia (Ridgway 2007), resulting in ocean warming that is among the top 10% of rates globally (Hobday and Pecl 2014).
Pelagic fishes were selected for analysis based on the availability of species occurrence records within the New South Wales Department of Primary Industries gamefish tagging database (NSW Department of Primary Industries 2019). This database contains geolocated (to 0.01 spatial resolution) and date-stamped (to a daily temporal resolution) occurrences of pelagic fishes tagged off eastern Australia by recreational anglers between 1974 and the present. Occurrence records were restricted to 1998-2018 to match the availability of satellite-derived oceanographic variables and filtered to ensure spatiotemporal independence among observations of the same species by only retaining data from a unique day and location as per Brodie et al. (2015) and Champion et al. (2018). Following these procedures, species with a minimum of 300 occurrences located nearshore of the continental shelf-break (i.e., 200-m isobath) were retained to enable robust statistical analysis. The resulting species analyzed were yellowtail kingfish (S. lalandi; n = 1088; hereafter "kingfish"), Australian bonito (Sarda australis; n = 314; hereafter "bonito"), Australian spotted mackerel (Scomberomorus munroi; n = 410; hereafter "spotted mackerel"), narrow-barred Spanish mackerel (Scomberomorus commerson; n = 889; hereafter "Spanish mackerel"), and common dolphinfish (Coryphaena hippurus; n = 558; hereafter "dolphinfish"). A sample of 300 occurrences were randomly selected from the total remaining number of species records for model selection and evaluation to ensure model results were comparable among species (Fig. 1A). Spatial and temporal variograms were generated to assess for spatiotemporal autocorrelation among the resulting set of species occurrence records (Fig. S1), which confirmed that the data processing methods applied resulted in spatially and temporally independent sets of species occurrence data. A sample of 10,000 pseudo-absences was randomly distributed nearshore of the continental shelf-break throughout the study extent and period to characterize unsuitable environmental habitat and generate a binomial response variable for statistical modeling (Barbet-Massin et al. 2012). A sensitivity analysis was undertaken to assess the relationship between the number of pseudo-absences used in model fitting and the predictive accuracy of SDMs. This analysis found that the predictive accuracy of SDMs peaked when between 8000 and 12,000 pseudo-absences were used in model fitting and did not increase when a greater number of pseudo-absences were used ( Fig. S2), therefore supporting the selection of 10,000 pseudoabsences for developing SDMs.

Oceanographic and topographic variables
Oceanographic variables known to influence the distributions of pelagic fishes (Hobday and Hartog 2014) were downloaded from the Copernicus Marine Environment Monitoring Service (https://marine.copernicus.eu) and matched to occurrence and pseudo-absence data. These were sea surface temperature, salinity, eddy kinetic energy, sea level anomaly, and chlorophyll a (Chl a) concentration (Table S1). To quantify the influence of topographic habitat associations on species distributions, occurrence and pseudo-absence data were also matched to an index of topographic heterogeneity derived from the General Bathymetric Chart of the Oceans dataset (GEBCO_2020; https://www.gebco.net). This index is defined as the coefficient of topographic variation (CTV; Fig. 1B; Amatulli et al. 2018), which is a unitless measure of topographic heterogeneity within a defined spatial extent, calculated as: where SD pixel is the standard deviation of groups of pixels bilinearly interpolated from a gridded resolution of $ 0.004 (i.e., the native resolution of the GEBCO_2020 dataset) to 0.1 (i.e., the common resolution of environmental covariates used in model spatial predictions), and Median extent is the median bathymetric value calculated nearshore of the continental shelf-break throughout the study extent (Costello et al. 2010). Our analysis was restricted to the continental shelf to (1) align with the spatial distribution of fishery-dependent species occurrence records and (2) avoid bathymetric data offshore of the continental shelf inflating the denominator (i.e., Median extent ) when calculating CTV and thus masking topographic variation within the continental shelf environment. This index of topographic heterogeneity was used to assess for species associations with areas of high topographic variability, which are indicative of the presence of topographic features within the seascape (Hobday and Campbell 2009).
Collinearity among environmental variables was assessed using Spearman rank correlation coefficients. When correlated (> 0.5 and < −0.5) variables were identified, the covariate with the clearest ecological interpretation from covarying pairs was retained for model selection (Zuur et al. 2013). Sea surface temperature and salinity were negatively correlated, therefore salinity was removed from the suite of potential oceanographic predictors prior to model fitting. Full descriptions of environmental covariates available for model selection are presented in Table S1.

Species distribution modeling
Generalized additive mixed effects models (GAMMs) were used to describe environmental habitat suitability for all species throughout their eastern Australian distributions. For each species, two SDMs were developed that (1) quantified species preferences for multiple oceanographic variables (hereafter referred to as "oceanography-only" models) and (2) incorporated species preferences for topographic heterogeneity in conjunction with preferences for multiple oceanographic variables (hereafter referred to as "topography-inclusive" models).
GAMMs used a logistic link function to relate binomial distributed species occurrence and pseudo-absence data to multiple environmental covariates. Oceanography-only models were designed to contain all significant oceanographic predictors of species occurrence from the full set of available covariates. Therefore, model selection for species' oceanography-only models was undertaken by including all oceanographic covariates in an exploratory model, dropping individual covariates that were not significant predictors of species occurrence (at alpha level = 0.05) in order of least significance and refitting the model to assess for the statistical significance of all remaining covariates. Oceanography-only models did not omit covariates that significantly affected the occurrences of study species and were therefore representative of the dynamic oceanographic variables that simultaneously drive the distributions of coastal-pelagic fishes of eastern Australia (Brodie et al. 2015). Species responses to topographic heterogeneity were then included within all optimal oceanography-only models, producing two SDMs for each species (Table 1;  Table S2). Calendar year was included as a random term in all models to account for unmeasured annual variation in species catch-per-unit-effort. Nonlinear responses were modeled with penalized regression spline smoothers using generalized crossvalidation to optimize smooth functions and to avoid overfitting to the data (Zuur et al. 2009). Smoothers were replaced with linear terms if their effective degrees of freedom were approximately equal to 1, which indicates linearity with log-of-odds transformed binomial responses variables.
Model accuracy and predictive skill was assessed by (1) applying a temporally explicit sixfold cross-validation approach that trained and tested models on alternative subsets of data from the periods 1998-2008 and 2009-2018 to estimate the mean area under the receiver operating characteristic curve (AUC) and mean true skill statistic (TSS) for all oceanography-only and topography-inclusive models, and (2) comparing day-and location-specific spatial predictions of environmental habitat suitability with a set of 50 independent occurrence records for each species. A full description of this model evaluation process is presented in the Supplementary Material.

Range shift analyses
Oceanography-only and topography-inclusive SDMs were used to spatially predict species environmental habitat suitability throughout the study extent at a monthly temporal resolution from January 1998 to December 2018. Spatial predictions that fell within the upper 20% of environmental habitat suitably were classified as species' core environmental habitat and determined the poleward edge of species' core distributions within the study extent through time. The selection of the upper 20% of habitat suitability values to represent species core environmental habitats aligns with previous research assessing range shifts in the core distributions for pelagic fishes off eastern Australia (Robinson et al. 2015). While the latitudinal location of the poleward edge of core environmental habitat is dependent on the selection of habitat suitability thresholds, comparisons of range shift rates remain robust provided that (1) all spatial predictions are from models trained on the same amount of species occurrence and (pseudo-) absence data and (2)  Rates of redistributions predicted by oceanography-only and topography-inclusive SDMs were quantified and compared using linear mixed effects models (i.e., range shift models). For individual species and for all species collectively, range shift models took the form (in script notation): where Latitude is the most poleward latitude within the study extent corresponding to the predicted distribution of species' core environmental habitats, modeled as a function of time (Year) and the presence or absence of topographic habitat associations within predicted species distributions (SDM type ). An interaction between Year and SDM type was included to test for differences in rates of species range shifts predicted by oceanography-only and topography-inclusive SDMs over the study period. The coefficient of the fixed effect of Year on the latitudinal location of species core environmental habitats was used to estimate range shift rates in kilometers per decade. Calendar month (Month) was included as a random intercept term to account for natural intra-annual climate variability on the latitudinal distribution of species core oceanographic habitats throughout the study period.
Statistical analyses were undertaken using the software R (R Core Team 2019), with GAMMs fitted using the "gamm4" package (Wood and Scheipl 2013), k-fold cross-validation undertaken using the "dismo" package (Hijmans et al. 2013) and linear mixed effects models fitted using the "lme4" package (Bates et al. 2015). Alpha level = 0.05 was used to denote statistical significance. Residual plots were assessed visually to confirm that all linear mixed effects models satisfied the assumptions of normality, independent errors and homogeneity of variance. Spatial data were plotted in MATLAB (ver. 9.8, The MathWorks).

Environmental habitat preferences
Topographic heterogeneity was a significant predictor of the distributions of all pelagic fishes investigated, except dolphinfish ( Fig. 2; Table S2). Topographic heterogeneity had a positive linear effect on the occurrence of kingfish and a positive, approximately exponential, effect on the occurrences of bonito and both spotted and Spanish mackerel (Fig. 2). Sea surface temperature and sea level anomaly were significant predictors of the occurrence of all study species, while Chl a concentration was a significant predictor of the occurrence of all species except bonito. Eddy kinetic energy was a Table 1. Summary of the optimal oceanography-only ( * ) and topography-inclusive (◎) SDMs for each study species. Smoothing factors applied to environmental covariates are indicated by "s." AUC and TSS model evaluation statistics are derived from temporally explicit sixfold cross-validation, where oceanographic-only and topographic-inclusive models for each species were trained and tested on subsets of data from the periods 1998-2008 and 2009-2018. SST, sea surface temperature; EKE, eddy kinetic energy; CHL, chlorophyll a concentration; SLA, sea level anomaly; TOPO, coefficient of topographic variation.

Champion and Coleman
Seascape topography slows marine range shifts significant predictor in bonito and dolphinfish SDMs but not spotted or Spanish mackerel models (Table S2).

Habitat model evaluation
Temporally explicit sixfold cross-validation revealed that oceanography-only and topography-inclusive SDMs for all species had fair to good predictive accuracy (i.e., AUC > 0.7 and TSS > 0.5; Table 1; Pearce and Ferrier 2000). Topographyinclusive SDMs were associated with greater AUC and TSS evaluation statistics than oceanography-only models for all species except dolphinfish, indicating that the predictive accuracy and skill of oceanography-only models was generally improved by the incorporation of species topographic habitat associations (Table 1). Mean environmental habitat suitability values matching the location of independent species occurrence exceeded 0.4 and 0.5 for oceanography-only and topography-inclusive SDMs, respective. Habitat suitability values corresponding to independent occurrence data were significantly greater when habitat predictions were created using topography-inclusive models relative to predictions created using oceanography-only models for all species except dolphinfish (Fig. S3).  1998-2007 and 2008-2018. Note that greater latitudinal shifts in the poleward edge of species core environmental habitat (denoted by red dashed lines) between time periods are produced by models that excluded species preferences for topographic habitat, relative to models that included these environmental habitat associations.

Effect of seascape topography on range shift rates
Linear mixed effects models revealed significant poleward shifts in the leading edge of core environmental habitat of all species, irrespective of whether topographic habitat associations were included in SDM spatial predictions ( Fig. 3; oceanography-only SDMs: int = 343.17, slope = −0.18, t = 7.67, p < 0.001; topography-inclusive SDMs: int = 199.13, slope = −0.11, t = 4.42, p < 0.001). However, including species' topographic habitat preferences in SDM predictions significantly reduced rates of poleward range shifts for all pelagic fishes analyzed ( Fig. 4F; difference in rates of predicted redistribution = 88.8 [AE 21.2 SE] km per decade; t = 2.07, p = 0.03). Rates of poleward redistributions in individual species were slowed by 5.3-60.4%, or reduced by up to 131.4 km per decade, by seascape topography. Reductions in rates of predicted range shifts due to the influence of topography were statistically significant for three out of five species and spatially meaningful for all pelagic fishes analyzed. Topographic habitat associations resulted in the greatest reduction in the predicted rate of poleward redistribution for bonito (reduced by 60.4% or 131.4 [AE 31.7 SE] km per decade; t = 2.08, p = 0.04; Fig. 4B Fig. 4D), and dolphinfish (reduced by 5.3% or 7.1 [AE 23.5 SE] km per decade; t = 0.10, p = 0.92; Fig. 4E). The minor effect of topographic habitat on rates of redistribution identified for dolphinfish relative to other species reflects that topographic heterogeneity was not a significant predictor of the distribution of this species (Fig. 2).

Discussion
Accurately predicting spatiotemporal changes in species distributions under global change relies on quantifying key relationships between species and their environment (McHenry et al. 2019;Brodie et al. 2020). Our results demonstrate that seascape topography can be a significant predictor of pelagic fish distributions, and incorporating these environmental responses into SDMs can improve model skill and accuracy relative to models that only contain dynamic oceanographic covariates. Incorporating the effect of seascape topography on pelagic fish distributions into range shift analyses significantly reduced predicted rates of climate-driven redistributions by an average of $ 89 km per decade relative to analyses based only on species preferences for dynamic oceanographic variables. This suggests that pelagic fishes tracking changes in preferred oceanographic conditions also , and (F) all study species' core environmental habitats (dashed lines denote 95% confidence intervals) using models that did (topography-inclusive models; black data) and did not (oceanography-only models; red data) incorporate species responses to topographic habitat. associate with static topographic features encountered throughout the seascape, and that dynamic and static environmental variables combine to influence rates of species range shifts. Range shifts analyses that omit species preferences for seascape topography are likely, therefore, to predict rates of climate-driven redistributions that are $ 60% more rapid than analyses that include topographic habitat associations.
While widely distributed pelagic fishes are typically considered inhabitants of open ocean and coastal-shelf environments (Sund et al. 1981), the spatial distributions of these species are not random but aggregated in response to environmental and biological drivers (Briscoe et al. 2016). Even among groups of fishes that are classified as having "pelagic" or "open water" distributions, seascape topography can have variable effects on species occurrence and local abundance. For example, Hobday and Campbell (2009) found that catch rates of kingfish, bonito and southern bluefin tuna (Thunnus maccoyii) were 26, 6.6, and 2.5 times higher, respectively, at topographic features than away from them, while skipjack tuna (Katsuwonus pelamis) catch rates were one-third lower near topographic features. These results demonstrate variable effects of topographic heterogeneity on species relative abundances, which are supported by our finding that the occurrence of four out of five study species responded positively to topographic heterogeneity, while the distribution of dolphinfish was not significantly affected.
Variation in the strength of species topographic habitat associations corresponded with the magnitude of reductions in predicted rates of redistributions. For example, rates of redistribution in the poleward edge of core environmental habitat for species with the strongest associations with topographic habitat (bonito and kingfish) were reduced by 60.4 and 58.3 km per decade, respectively, while a reduction of 7.1 km per decade was identified for dolphinfish, which had a weak association with seascape topography. These findings are supported by the positive effect of topographic features on the relative abundance of bonito and kingfish found by Hobday and Campbell (2009) and Rees et al. (2018), and suggest that pelagic fishes that do not exhibit these habitat preferences may be undergoing more rapid range shifts than fishes that actively associate with topographic features. Thus, the effect of topographic habitat on pelagic fish distributions may be an important driver of variation in rates of climatedriven range shifts, which may manifest in punctuated rather than gradual shifts in species distributions depending on the presence or absence of topographic features throughout the seascape. Importantly, our results are based on rates of change in the poleward, or leading, edge of species' core environmental habitats, and it remains uncertain whether the centroid or trailing edge of species distributions are similarly affected by topographic heterogeneity.
Range shift analyses for pelagic fishes commonly focus on dynamic oceanographic habitats due to the availability of satellite-measured oceanographic predictors of pelagic species distributions (Hobday and Hartog 2014;Champion et al. 2019). We demonstrate how a user-specified index of topographic heterogeneity can be used to incorporate the effect of the physical seascape into predictions of species redistributions. Bathymetric data are freely available and have global coverage, allowing for species associations with seascape topography to be quantified in any region of the ocean (Amatulli et al. 2018;Costello et al. 2010). Improved efforts to capture important yet often overlooked species-environment relationships, for example, as additional covariates in correlative SDMs as we have done, are crucial given their significant effect on predictions of species distributions (McHenry et al. 2019). Doing so will improve estimates of the rate at which climate change is altering the distribution of life on earth, which may be biased if species associations with static environmental variables are not incorporated into range shift analyses.