Global habitat loss of a highly migratory predator, the blue marlin (Makaira nigricans)

Climate change is driving the redistribution of species throughout the oceans. However, the speed and magnitude of species responses, including shifts in their distribution, are variable and species specific. Quantifying the effect of environmental conditions on species distributions is crucial to informing management and conservation efforts. Blue marlin (Makaira nigricans) is a wide‐ranging top predator occurring circumglobally in tropical and subtropical waters and is heavily impacted by international longline fisheries. This study aimed to predict the global distribution of blue marlin and characterize the effects of climate variability thereon.


| INTRODUC TI ON
Understanding species habitat preferences and how these preferences influence their distributions and habitat use is increasingly important due to the accelerating effects of global climate change (Burrows et al., 2011;Pecl et al., 2014;Poloczanska et al., 2013).
Shifts in species distributions are occurring substantially faster in marine systems compared to terrestrial (Chen et al., 2011;Poloczanska et al., 2013;Sorte et al., 2010) and are primarily shifting in poleward directions (Champion et al., 2018;Erauskin-Extramiana et al., 2019;Hill et al., 2016;Sorte et al., 2010). Responses to climate change, however, vary substantially among species (Pinsky et al., 2013) and are strongly influenced by species-specific traits (Champion et al., 2021;Sunday et al., 2015). In oceanic systems, highly mobile species may have large annual migrations brought about by their ability to follow preferential environmental conditions Fredston et al., 2021;Sunday et al., 2015). Therefore, it is important to better understand the habitat preferences of highly mobile species. This information provides a more informed understanding of the drivers of their distribution and greater insight into the measures needed to manage and conserve species facing climate-induced range shifts.
Pelagic fish species such as tunas and billfishes are highly mobile and often undertake long distance migrations. Tagging studies have shown that billfishes and tunas have annual migrations that can traverse ocean basins and tens of thousands of kilometres to reach spawning and/or foraging grounds (Block et al., 2005;Ortiz et al., 2003). These species are targeted and caught as bycatch in numerous longline fisheries and recreational fisheries around the world (Arocha & Ortiz, 2006;FAO, 2020;Graves & Horodysky, 2010;Ortiz & Farber, 2001). As these fish are epipelagic, upper trophic level predators, they also play important roles in maintaining the health and function of marine ecosystems (Myers et al., 2007;Myers & Worm, 2003;Scheffer et al., 2005). Quantitatively characterizing the distributions and habitat use of these top predators and understanding their potential responses to climate change is therefore vital for timely and effective management and conservation decisions.
Blue marlin (Makaira nigricans) is a large (maximum weight > 700 kg; Hill et al., 1989), highly migratory species of billfish (family Istiophoridae) distributed globally in tropical and subtropical waters. Although known to spawn in several locations such as Japan, Hawai'i, French Polynesia, the Straits of Florida and Bahamas (Hopper, 1990;Howard & Ueyanagi, 1965;Luthy et al., 2005;Richardson et al., 2009;Serafy et al., 2004;Shimose et al., 2009Shimose et al., , 2012, the drivers of blue marlin movements are still poorly understood. Blue marlin is frequently caught as bycatch in tuna longline fisheries (Molony, 2005;Serafy et al., 2004) but also targeted in artisanal and recreational fisheries (Arocha & Ortiz, 2006;Brinson et al., 2006;Luckhurst, 2003). The combined catches of these fisheries have led to a listing of vulnerable by the International Union for Conservation and Nature (Collette et al., 2011). Blue marlin is considered a single global species (Collette et al., 2006) with separate stocks in the Atlantic, Pacific and Indian Oceans (Graves & McDowell, 2001;Hinton, 2001;Williams et al., 2020), although information on blue marlin in the Indian Ocean is limited (Graves & McDowell, 2015). However, recent genetic evidence suggests blue marlin in the Pacific and eastern Indian Oceans constitute a single stock (Chen et al., 2016).
Acoustic and satellite studies indicate that blue marlin is epipelagic with a bimodal depth distribution (Carlisle et al., 2017;Goodyear et al., 2008). Daytime habitat is primarily between 25 and 100 m with periodic excursions into deeper water, whereas nighttime habitat is generally within the top 10 m of the water column (Block et al., 1992;Carlisle et al., 2017;Goodyear et al., 2008;Holland et al., 1990;Kraus et al., 2011). Sea surface temperature (SST) and dissolved oxygen are thought to be significant drivers of habitat use for blue marlin. This species is generally found at temperatures above 24°C but it prefers warmer waters up to 30°C (Block et al., 1992;Carlisle et al., 2017;Goodyear et al., 2008;Graves et al., 2002;Holland et al., 1990).
Shoaling of the oxygen minimum zone in the Eastern Tropical Pacific (ETP) and Atlantic (ETA) are thought to limit the vertical habitat of blue marlin (Prince et al., 2010;Prince & Goodyear, 2006). Low temperature and low dissolved oxygen have a similar effect when they co-occur in the Pacific (Carlisle et al., 2017). Environmental variability related to the El Niño Southern Oscillation (ENSO) cycle is also thought to influence the habitat use of blue marlin in the Pacific. Su et al. (2011) found that the distribution of blue marlin shifted eastwards during the 1997-1998 El Niño associated with increasing SSTs and a deepening of the thermocline. Carlisle et al. (2017) suggested the westward extension of cold, low oxygen water into the Central Pacific associated with the 2010 La Niña may act as a barrier to transequatorial migrations into the South Pacific. Haulsee et al. (2022) found negative correlations between blue marlin sightings from recreational logbook data and ENSO events in more localized areas of the ETP. Together, these studies suggest that blue marlin distribution is strongly influenced by the physical environment.
Species distribution models (SDMs) are important tools that are increasingly being used to quantify the spatial distribution of species habitat in terms of environmental conditions (Block et al., 2011;Elith & Leathwick, 2009;Guisan & Thuiller, 2005;Robinson et al., 2017).
Previous studies have used SDMs to evaluate the effect of environmental variables and seasonal variation on blue marlin distribution and abundance in the Pacific using longline fisheries data (Su et al., 2008(Su et al., , 2011. The distribution of blue marlin in the Atlantic was inferred using temperature and depth data from electronic tagging (Goodyear, 2016). However, these studies did not incorporate tagged individuals' horizontal movement data and were specific to K E Y W O R D S blue marlin, global, habitat, loss, migratory, predator, shift, species distribution model a single ocean basin. Despite the continued fishing pressure and potential long-term impacts of climate change on habitat availability, an evaluation of their potential habitat distribution is lacking at the global scale. In this study, the habitat preference of blue marlin is evaluated, and its distribution is predicted from 2000 to 2016 throughout their worldwide range using a unique, long-term telemetry data set, collected mostly through the International Game Fish Association Great Marlin Race (IGMR). Seasonal variability of blue marlin habitat and the potential effects that climate variability may have had on its distribution is also investigated.

| Species data
Data on the movements of blue marlin were obtained using popup archival transmitting (PAT) tags, (MK10 and miniPAT tags, Wildlife Computers) deployed around the world between August 2000 and February 2017 ( Figure 1). The programming scheme for PAT tags changed over the years as tag technology and analytical approaches improved, with tags being programmed to release initially after 120 days and subsequently for 180, or 240 days. Critical to the deployments was development of a new tether system that helped to keep the tags attached for longer durations than earlier studies. This system consisted of a custom-built titanium dart, a 3layer tether constructed of monofilament, a sleeve of Kevlar line that increased abrasion resistance and 1/8th inch shrink wrap (Wilson et al., 2015). Stainless steel and titanium pole application tips were set to 15-18 cm lengths and built with deep set insets to help create a system for deeper application into the blue marlin hypaxial muscle upon initial tagging. Tags recorded pressure, temperature and light data. Daily geolocation estimates were generated from light levels, SSTs recorded from the tags, and remotely sensed SST following Teo et al. (2004). Light and SST-based geolocation tracks were then refined using a Bayesian state-space model (SSM), which also estimated the uncertainty around daily location estimates (Block et al., 2011;Winship et al., 2012). Median values of tag deployment data in the dataset (i.e., deployment duration, distance travelled, travel rate and fish weight) were tested for differences both among (Atlantic, Pacific, Indian) and within (northern vs. southern hemisphere) oceans using a Kruskal-Wallis One-Way Analysis of

| Environmental data
Ten environmental variables were selected to describe the habitat preference of blue marlin based on their importance in previous studies (Brodie et al., 2018;Carlisle et al., 2017;Madigan et al., 2021): SST, surface chlorophyll-a (Chl-a), sea surface height (SSH), oxygen at 100 m depth, total kinetic energy (TKE), mixed layer depth (MLD), bathymetry, rugosity and the standard deviation of SST and SSH (an index of frontal activity, Brodie et al., 2018) (Table S1). Oxygen at 100 m was chosen as a 2-dimensional representation of oxygen depth limitation (Carlisle et al., 2017), which may influence the horizontal distribution of blue marlin. Daily environmental data for tag locations were extracted, with data averaged over a spatial area that matched the size of tag error estimates ( Figure S1). Tag data with error estimates larger than 3° were removed from analysis. Tag F I G U R E 1 Map of state-space-modelled blue marlin tracks used in the study (n = 144). Tracks are colour coded by year and range from 2000-2017 error estimates of 0 were averaged at a fixed spatial resolution of 0.25° (Hazen et al., 2016). For standard deviation variables (SST_sd, SSH_sd), the same approach was used but the standard deviation over the spatial area defined by the tag error estimates was calculated. Correlations between environmental variables were examined for strong relationships using the Spearman Rank method prior to model development. Values for Chl-a, TKE, MLD, SST_sd, SSH_sd and rugosity were log-transformed prior to inclusion in modelling due to the presence of excessive skew in their distributions.

| Species distribution model
A generalized additive mixed model (GAMM) was used to estimate global blue marlin habitat suitability. Presence/pseudo-absence was modelled as a function of environmental variables using a binomial family, with spline-smoothers fit to each environmental covariate and individual tag id as a random effect. Each tag daily location estimated using SSM was used as presence data in the model.
Pseudo-absences were simulated using Correlated Random Walks (CRWs) and were used to allow binomial modelling of species habitat suitability (Hazen et al., 2016. For each track, 100 CRWs of equal length and starting from the tagging location were generated using randomly selected step lengths and turning angles of each tag (Hazen et al., 2016. The influence of the random pseudo-absence CRW choice on model selection was then assessed (Hazen et al., 2016). The selected model was run 60 times on the full dataset, each with a unique set of randomly selected CRW pseudo-absence tracks, and the percentage of the 60 model iterations each covariate was significant was calculated. K-fold cross validation was used to determine the best CRW set to use with the final model. The full dataset was divided into five folds, with four of the subsets combined to train the model and the fifth used as a test subset. The selected model was then run 60 times on the training subset, each with a unique set of randomly selected CRW pseudo-absence tracks and the AUC from model predictions on the test subset was calculated. This process was repeated five times so that each subset was used as a test set. The CRW set with the highest mean AUC was selected for use in the final dataset. GAMMs were built and evaluated using the MGCV (version 1.8-28) and ROCR (version 1.0-7) packages in the R Statistical Environment (version 4.1.2, R Core Team, 2020). Marginal means for each covariate were calculated using the R package ggeffects (version 1.1.2; Lüdecke, 2018). The final model with the selected covariates was then fit to the full dataset which included the best performing set of CRWs.
Fisheries-dependent catch-per-unit-effort (CPUE) data were used as an additional validation of model performance. CPUE data were obtained from four Regional Fishery  Nominal CPUE was calculated as the ratio of blue marlin catch (tons) to the number of hooks for longline data aggregated from 2000 to 2016 at 5° × 5° spatial resolution. The correlation between predicted habitat values (averaged to 5° × 5° spatial resolution) and CPUE was calculated using a modified t-test from the R package "SpatialPack" (version 0.3-8196; Osorio et al., 2016), which assumes correlations to be inflated by spatial autocorrelation.
The best-fit model consisting of 10 environmental variables and fit with the best-performing set of CRW pseudo-absences was used to predict daily global blue marlin habitat from 2000 to 2016 at a spatial resolution of 0.25° × 0.25°. Daily predictions were then averaged to a monthly scale. Although model fitting included data from 2017, there were only a few tracks of short duration during this year. Therefore, habitat predictions were not made beyond 2016.
In this framework, a single model fit to the entire dataset was used to predict seasonal trends in habitat suitability rather than fitting models to specific seasons or months. The degree to which predictions were extrapolated beyond the environmental space used for model fitting was assessed. Extrapolation was quantified using the compute_extrapolation function in the dsmextra R package (version 1.1.5; Bouchet et al., 2020). For each month from 2000 to 2016, the percent of predicted grid cells characterized by environmental conditions similar to those experienced in the model fitting dataset was calculated. Higher percentages indicate less extrapolation.

| Spatial and time-series analyses
Trends in habitat suitability were evaluated based on monthly median habitat values. For each month, the predicted habitat space was subset based on the optimal cutoff calculated with ROCR into high suitability (predicted habitat > 0.52) and low suitability (predicted habitat < 0.52). The median suitability was calculated for each subset, creating two time-series of monthly median habitat values from 2000 to 2016. Yue & Wang (2002) time-series methods were used to test for a significant trend in median values using the R package "spa-tialEco" (version 1.3-7; Evans, 2020). The test generates a Kendall's τ statistic ranging from −1 to 1 with values closer to 0 indicating no significant trend. Additionally, the monthly habitat predictions were divided into two equal groups representing the first and second halves of the study period from 2000 to 2008 and 2009 to 2016.
The first half was subtracted from the second half to identify the spatial scope of habitat change over the course of the study.
Core habitat of blue marlin and estimated spatial anomalies in their habitat suitability over time were identified. Core habitat was defined as the top 25th percentile of habitat values calculated from the mean habitat of the study period represented as the mean of all months across all years. Interannual variation was estimated from yearly mean habitats represented as the mean of all months within a given year. Intermonth variation was estimated from monthly mean habitats represented as the predicted habitat of a given month averaged over all years. Monthly habitat anomalies were calculated as the monthly mean habitat subtracted from the mean habitat of the study period. Negative anomalies represent below average values and positive anomalies represent above average values.
Declines were most evident in the North Atlantic, the majority of the Indian Ocean, the South Pacific east of French Polynesia

| DISCUSS ION
Blue marlin is an epipelagic scombroid fish species that is biomechanically designed for long distance travel (Long, 1989). Our understanding of blue marlin distribution and habitat preference has been limited by the difficulty of collecting fishery-independent data due to their oceanic distribution and highly migratory behaviour F I G U R E 2 GAMM response curves for blue marlin probability of presence model. Probability of presence on the y-axis calculated from the marginal means of each variable. Predictor variable range on the x-axis (actual data shown as rug plots). Shaded regions are the 95% confidence intervals. Chla, chlorophyll-a; MLD, mixed layer depth; SSH, sea surface height; SSH_sd, sea surface height standard deviation; SST, sea surface temperature; SST_sd, sea surface temperature standard deviation; TKE, total kinetic energy. Log values are calculated as natural logarithms (Collette et al., 2011). The development and advancement of satellite tag technology in the past two decades has greatly enhanced the collection of fishery-independent data across a broad range of marine taxa (e.g., Block et al., 2011), however billfishes have remained challenging in terms of tag retention due to premature detachment issues. Additionally, the expense to purchase and deploy tags requires significant financial and time investment (e.g., Lutcavage et al., 2015). This study represents a highly successful application of citizen-based science through a collaboration between Stanford demonstrate the large-scale decline in suitability of core habitat as well as the concurrent increase in suitability of marginal habitat at the latitudinal edges of the range. The tradeoff of this approach is that some of the dynamics within specific regions of the different ocean basins will necessarily be lost once integrated into global patterns. Some of these dynamics may play important roles regionally, and warrant further investigation, but that is beyond the scope of this paper.
Model validation with an independent dataset is considered the gold standard in species distribution modelling (Araújo et al., 2019).
Although direct comparison of model predictions and relative CPUE data has limited value due to the inherent differences between telemetry and fishery-dependent data (Brill & Lutcavage, 2001;Campbell, 2015;Maunder et al., 2006;Maunder & Punt, 2004), a positive correlation provided additional validation of the model in this study. Consistent with previous smaller studies on this species, predicted suitable habitat for blue marlin was centred on the equator and spanned a broad latitudinal range (Goodyear, 2016;Marín-Enríquez et al., 2020;Su et al., 2008). Over the course of the study, there was little variation in predicted habitat between years despite an overall decline in habitat suitability. In the Pacific, Su et al. (2011) found Likelihood of habitat and standard deviation ranges from low (white and blue) to high (orange and red) with strong associations above 23°C, consistent with temperature preferences found in previous studies (Block et al., 1992;Carlisle et al., 2017;Goodyear et al., 2008;Graves et al., 2002;Holland et al., 1990;Su et al., 2011). In the first half of the year (January to July, austral summer and fall), habitat suitability at the lower latitudes near the edge of the southern habitat were above average, whereas habitat near the edges of the norther latitudes were below average.
This pattern becomes more neutral in July and fully reverses through the second half of the year (July to December, boreal summer and fall). A similar pattern derived from catch-effort data was found in the Atlantic (Goodyear, 1999(Goodyear, , 2003, with a northward/southward expansion of CPUE in the summer months. Blue marlin in the Pacific was predicted to move north from April to October when predicted habitat suitability increases and south thereafter (Su et al., 2011).
That pattern coincides with observed North to South equatorial crossings of PAT tagged fish from Hawai`i occurring between late September and mid-November (Carlisle et al., 2017) and the decline in habitat suitability in the South Pacific from July to December.
Oxygen at 100 m was the third most important variable in the model, however its impact on predicted habitat suitability appears to be primarily at the latitudinal extremes of the habitat (Figure 3).
Positive associations in blue marlin presence were found where oxygen at 100 m ranged broadly from 0.54 and 5.1 ml L −1 . Values below 3.5 ml L −1 are generally considered to be stress inducing (Bushnell & Brill, 1991;Ingham et al., 1977;Prince & Goodyear, 2006) and values at the low end of the range considered lethal (Brill, 1994). However, predicted habitat suitability remained elevated in areas where oxygen at 100 m fell below these "stressful" levels, such as the ETP Mature blue marlin are thought to migrate between foraging and spawning habitats (Shimose et al., 2009;Shimose et al., 2012).
Limited vertical habitat could have less impact on horizontal space use if occurring during more directed, migratory movements where SST may be the driving environmental factor (e.g., Carlisle et al., 2017;Su et al., 2011). However, regions where low dissolved oxygen could limit vertical habitat may also provide increased foraging opportunities. Prince and Goodyear (2006) suggested these regions could lead to enhanced interactions between blue marlin and its prey due to aggregation of prey into the same narrow mixed layer. They also found a sympatric species, sailfish (Istiophorus Highly suitable habitat, including 96% of core habitat, declined in suitability throughout most of the range in all three oceans over the course of the study. During the same period, the suitability of marginal habitat increased, primarily at the edges of the range ( Figure 6).  Triacca, 2021). Several studies have noted habitat shifts and predicted poleward expansion of fish habitat based on climate change models through the end of the century (Champion et al., 2018(Champion et al., , 2021Erauskin-Extramiana et al., 2019;Fredston-Hermann et al., 2020;Hazen et al., 2013). In this study, the present-day loss of highly suitable habitat throughout much of the range with a concurrent poleward expansion of marginal habitat is demonstrated, similar to what has been reported for black marlin in the southwest Pacific (Hill et al., 2016). Concerningly, this result suggests that ocean warming due to climate change may be making equatorial waters less suitable even to highly mobile, tropical/subtropical species.
As a highly migratory species, blue marlin is likely to adapt to environmental changes by following preferred habitat as it shifts, especially if the primary prey are experiencing similar shifts in habitat (e.g., Carroll et al., 2019;Erauskin-Extramiana et al., 2019). Using climate models to project blue marlin habitat into the future could help provide insight into the extent of these range shifts and inform the management and conservation policies needed to ensure the longterm viability of this species. Importantly, mapping highly migratory fish movements now and in the future will help to assess interactions with fisheries fleets setting longlines for more temperate species and anticipate interactions with blue marlin as it ranges increasingly into these habitats. Targeted and dynamic management policies such as time-area closures (Maxwell et al., 2020) and gear restrictions (e.g., Graves et al., 2012) may be most effective in areas where high levels of overlap are identified.

ACK N OWLED G EM ENTS
The Foundation and Stanford University.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Tracking data used in this study are available in Dryad data set https://doi.org/10.5061/dryad.rfj6q 57c6.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13606.