Ensemble Modelling of Skipjack Tuna (Katsuwonus pelamis) Habitats in the Western North Pacific Using Satellite Remotely Sensed Data; a Comparative Analysis Using Machine-Learning Models

To examine skipjack tuna’s habitat utilization in the western North Pacific (WNP) we used an ensemble modelling approach, which applied a fisherderived presence-only dataset and three satellite remote-sensing predictor variables. The skipjack tuna data were compiled from daily point fishing data into monthly composites and re-gridded into a quarter degree resolution to match the environmental predictor variables, the sea surface temperature (SST), sea surface chlorophyll-a (SSC) and sea surface height anomalies (SSHA), which were also processed at quarter degree spatial resolution. Using the sdm package operated in RStudio software, we constructed habitat models over a 9-month period, from March to November 2004, using 17 algorithms, with a 70:30 split of training and test data, with bootstrapping and 10 runs as parameter settings for our models. Model performance evaluation was conducted using the area under the curve (AUC) of the receiver operating characteristic (ROC), the point biserial correlation coefficient (COR), the true skill statistic (TSS) and Cohen’s kappa (k) metrics. We analyzed the response curves for each predictor variable per algorithm, the variable importance information and the ROC plots. Ensemble predictions of habitats were weighted with the TSS metric. Model performance varied across various algorithms, with the Support Vector Machines (SVM), Boosted Regression Trees (BRT), Random Forests (RF), Multivariate Adaptive Regression Splines (MARS), Generalized Additive Models (GAM), Classification and Regression Trees (CART), Multi-Layer Perceptron (MLP), Recursive Partitioning and Regression Trees (RPART), and Maximum Entropy (MAXENT), showing consistently high performance than other algorithms, while the Flexible Discriminant Analysis (FDA), Mixture Discriminant Analysis (MDA), Bioclim (BIOC), Domain (DOM), Maxlike (MAXL), Mahalanobis Distance (MAHA) and Radial Basis Function (RBF) had lower performance. We found inter-algorithm variations in predictor variable responses. We conclude that the multi-algorithm modelling approach enabled us to assess the variability in algorithm performance, hence a data driven basis for building the ensemble model. Given the inter-algorithm variations observed, the ensemble prediction maps indicated a better habitat utilization map of skipjack tuna than would have been achieved by a single algorithm.


Introduction
Species distribution modelling is based on the ecological niche concept, which is considered as the volume in the environmental space that permits positive growth of a species [1,2]. The ecological Remote Sens. 2020, 12, 2591 3 of 15 future conditions [26,27]. From this background, we found it feasible to apply a multi-algorithm modelling approach to derive skipjack tuna habitats in the western North Pacific, where to the best of our knowledge such an approach has not been employed for skipjack tuna with the number of algorithms we have used. Our choice of the sdm package was informed by its ability to execute a suite of 17 algorithms using RStudio, and allow for substantial flexibility in evaluation of model performance, predictions and generation of ensemble models [9,15].
The skipjack tuna (Katsuwonus pelamis) fishery in the western North Pacific is closely associated with the seasonal migration of the fish northwards from spring to summer, and southwards at the onset of winter, which is closely linked to feeding [28,29]. The fish track highly productive areas around oceanographic features observable from sea surface temperature (SST) and ocean color gradients, eddies and warm streamers [8,30,31]. They inhabit the upper mixed layer [32,33], and are opportunistic predators which feed voraciously on pelagic fishes, squids, a variety of crustaceans and young skipjacks [34,35]. They are known to migrate over large areas in search of high concentrations of forage [32]. Consequently, these bio-physical oceanographic signals with strong trophic linkages are important indicators of skipjack tuna aggregation sites. The signals are easily monitored with satellite observations and are, therefore, very useful for modeling skipjack tuna habitat formations through correlative or machine learning techniques. Various studies have provided detailed information on skipjack tuna habitat, both from surface [5,8,36] and sub-surface variables [37,38]. The role of temperature, chlorophyll-a concentration and currents in skipjack tuna habitat characterization is well known from several studies [21,29]. For instance, a lower limit SST at 18 • C, surface chlorophyll-a at 0.2-0.3mg/L and warm core eddies in the western North Pacific have been shown to indicate frontal oceanographic features that define skipjack tuna habitat [8,30,36].
Much of the work done on skipjack tuna habitat in the western North Pacific either using surface satellite data or sub-surface data relies heavily on single algorithm computations to relate tuna occurrence to their environment [5,29,30,36,37,39]. These studies are fundamental in elucidating skipjack tuna's habitat utilization, and can be enriched by exploring algorithms which have received little attention in this field. While multi-algorithm ensemble habitat predictions have been tried for other species in the North Pacific [22], we did not find work that has explored skipjack tuna habitat utilization in the western North Pacific using multi-algorithm species distribution models, hence our work seeks to enhance previous studies in this area [5,8,36]. Our specific objectives were: (i) to model skipjack tuna habitats in the western North Pacific using a suite of 17 presence-only species distribution algorithms and remotely sensed information; (ii) compare the performances of these algorithms using model performance metrics; and (iii) generate ensemble predictions of habitat suitability maps using the sdm package.

Study Area
The skipjack tuna caught in the western North Pacific (WNP) off Japan (18-50 • N and 125-180 • E) arrive in the area via the Nansei Islands and the Pacific Coast of Japan [40]. The oceanographic features in this area influence the migration routes ( Figure 1) of the fish as they migrate north [29,41]. The Kuroshio Current, which assumes three major paths south of Japan, is of critical importance to the skipjack's migration because it transports warm oligotrophic waters, forming pathways for migration, hence influencing formation of pelagic fisheries [42]. The behavior of the Kuroshio Extension, warm streamers and warm core rings (WCR) in the productive Kuroshio Oyashio Transition Zone (KOTZ) is an important feeding and migratory area for skipjacks, which also influences fishing ground formation [30,31]. Further north, towards the northern limit of skipjack tuna migration, the cold Oyashio Current waters flow southward [43], transporting low temperature, low salinity and nutrient rich waters to the sub-tropical gyre [44], and forming two southward tongue-shaped intrusions off Honshu, known as the First and Second Oyashio Intrusions [42,45]. These cold Oyashio waters and the Oyashio Front [46] inhibit skipjack tuna's northern migration because they avoid waters whose SSTs are below 18 • C [29]. The Oyashio meanders are separated by a WCR originating from the northward movement of the ring produced by the Kuroshio [47], and the frontal zones between the WCR and the meanders are key foraging areas [8]. The behavior of the major currents also influences the lower trophic level dynamics through the distribution and circulation of nutrients, hence the density of phytoplankton and zooplankton populations in the WNP [42]. The high densities of phytoplankton in high nutrient waters subsequently support large populations of various zooplankton species, which are fed upon by smaller nekton [43,46,48]. This creates a downstream trophic chain where skipjack tuna and other pelagic predators exploit such areas of high productivity to forage on the small organisms (squids, crustaceans, and fishes) [7,29,31,32]. Consequently, fishers of skipjack tuna in the WNP often track oceanographic features such as thermal and ocean color fronts, upwelling zones, and edges of large eddies [5,49], to locate areas with dense aggregations of skipjack tuna.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 17 because they avoid waters whose SSTs are below 18 o C [29]. The Oyashio meanders are separated by a WCR originating from the northward movement of the ring produced by the Kuroshio [47], and the frontal zones between the WCR and the meanders are key foraging areas [8]. The behavior of the major currents also influences the lower trophic level dynamics through the distribution and circulation of nutrients, hence the density of phytoplankton and zooplankton populations in the WNP [42]. The high densities of phytoplankton in high nutrient waters subsequently support large populations of various zooplankton species, which are fed upon by smaller nekton [43], [46], [48].

Fishery Data
Fishery data are often used as indicators of a species presence. We obtained daily skipjack tuna catch data from the Ibaraki Prefecture Fisheries Research Station, covering March to November 2004, the period when the skipjack tuna fishery is active in this area. These data were recorded from pole and line vessels fishing for skipjack tuna in the study area, which usually record their fishing locations and catches. The data were digitized, compiled into monthly composites and converted into 0.25 o resolution grids [50] in cells where catches were positive. Re-gridding the data was necessary to ensure that the fishery data matched the resolution of the environment grids. Matching the gridded presence data to predictor variables for a corresponding time period provides a record of the environmental conditions prevailing at the location where the fish were caught.

Environmental Data
A monthly predictor variable database was compiled consisting of SST, sea surface chlorophylla (SSC), and sea surface height anomalies (SSHA). Given their importance as environmental predictors of skipjack tuna habitat [28], [41], the three variables are widely used in modelling studies. Skipjack tuna have a wide tolerance range for ambient temperature, hence SST data are an important indicator of their distribution patterns [32], [41]. The SSC data provide information on ocean productivity [51], and are important for detecting fronts and eddies that are not always evident in SST maps [52], [53]. The elevated productivity around these features attracts large schools of tuna, which aggregate around them to feed on lower trophic level organisms [6], [31]. The SSHA data are an indicator of ocean dynamic topography, which provides information on movement of water

Fishery Data
Fishery data are often used as indicators of a species presence. We obtained daily skipjack tuna catch data from the Ibaraki Prefecture Fisheries Research Station, covering March to November 2004, the period when the skipjack tuna fishery is active in this area. These data were recorded from pole and line vessels fishing for skipjack tuna in the study area, which usually record their fishing locations and catches. The data were digitized, compiled into monthly composites and converted into 0.25 • resolution grids [50] in cells where catches were positive. Re-gridding the data was necessary to ensure that the fishery data matched the resolution of the environment grids. Matching the gridded presence data to predictor variables for a corresponding time period provides a record of the environmental conditions prevailing at the location where the fish were caught.

Environmental Data
A monthly predictor variable database was compiled consisting of SST, sea surface chlorophyll-a (SSC), and sea surface height anomalies (SSHA). Given their importance as environmental predictors of skipjack tuna habitat [28,41], the three variables are widely used in modelling studies. Skipjack tuna have a wide tolerance range for ambient temperature, hence SST data are an important indicator of their distribution patterns [32,41]. The SSC data provide information on ocean productivity [51], and are important for detecting fronts and eddies that are not always evident in SST maps [52,53].
The elevated productivity around these features attracts large schools of tuna, which aggregate around them to feed on lower trophic level organisms [6,31]. The SSHA data are an indicator of ocean dynamic topography, which provides information on movement of water masses, and by extension the flow of heat and nutrients, which subsequently influence productivity [54]. We downloaded daily resolved optimally interpolated SST global data from the National Oceanic and Atmospheric Administration's (NOAA) National Climatic Data Center (NCDC) ( Table 1), for the year 2004. The optimally interpolated data from the Advanced Very High Resolution Radiometer (AVHRR) and Advanced Microwave Scanning Radiometer on the Earth Observing System (AMSR-E) provided better coverage of the fishing locations because the effect of missing data due to clouds is eliminated. Monthly Aqua-MODIS (Moderate Resolution Imaging Spectroradiometer)~4 km standard mapped images of SSC for 2004 were downloaded from the National Aeronautics and Space Administration (NASA) Ocean Color data portal ( Table 1). The monthly averaged SSHA data were downloaded from the BloomWatch 180 portal ( Table 1). Data processing (monthly averaging and sub-setting) and mapping was undertaken in the SeaWiFS (Sea-viewing Wide Field-of-view Sensor) Data Analysis System (SeaDAS) [55]

Habitat Modeling Using sdm Package
The sdm package is based on a framework that uses presence-only data to compute habitat suitability indices (HSI), where the values closer to or equal to 1 represent the high potential habitat areas, and those close to zero represent poor habitat areas. The package comprises a suite of 17 algorithms, whose operation in RStudio provides an ideal platform with sufficient flexibility to adjust for algorithm selection, the associated parameter settings and generation of ensembles. We constructed 9 monthly habitat models in RStudio (version 3.6.1) using the 17 algorithms implemented in the sdm package [15]. We created a script in RStudio which loads all the necessary R packages as described by [15], ingests the fishery data and environmental variable ASCII grids, and subsequently calls the respective algorithms for model fitting, performance evaluation, algorithm selection and ensemble predictions. The skipjack tuna presence data were used as the presence variable and the SST, SSC and SSHA as predictor environmental variables. We applied a 70:30 split for training and test data respectively [57,58], with replication by bootstrapping, for 10 runs. In the bootstrapping technique, sampling with replacement is repeated, each time a sample with equal size as the original data is drawn and used for training data. The observations that are not selected are used for the evaluation at each run [15]. We used variable importance plots to assess the contribution of each predictor variable to model fit, and the response curves to show the relationship between the probability of occurrence of skipjack tuna and each of the predictor variables. Table 2 provides a summary of the 17 algorithms used and their dependent packages in R.

Evaluation of Model Performance
We evaluated the performance of our models using a multi-metric approach consisting of the area under the curve (AUC) of the receiver operating characteristic (ROC), the point biserial correlation coefficient (COR), the true skill statistic (TSS) and Cohen's kappa (k) metrics, which is an approach that has been applied in other studies [22,59]. The AUC metric measures the ability of a model to discriminate between sites where a species is present, versus those where it is absent, which provides an indication of the usefulness of the models for prioritizing areas in terms of their relative importance as habitat for the particular species [60]. The AUC ranges from 0 to 1, where a score of 1 indicates perfect discrimination, a score of 0.5 implies predictive discrimination that is no better than a random guess, and values <0.5 indicate performance worse than random [60,61]. The COR is the correlation between the observation in the presence-absence (pseudo-absences) dataset (a dichotomous variable) and the prediction, and is calculated as a Pearson correlation coefficient [61]. It is similar to AUC, but also accounts for how far the prediction varies from the observation, which gives further insight into the distribution of the predictions [61]. The TSS, calculated as sensitivity plus specificity less one, is an improvement of the k which takes into account sensitivity and specificity, and is insensitive to prevalence [62]. Like kappa, it takes into account both omission and commission errors, and success as a result of random guessing, and ranges from −1 to +1, where +1 indicates perfect agreement and values of zero or less indicate a performance no better than random. [62]. The Cohen's kappa is a widely used statistic which corrects the overall accuracy of model predictions by the accuracy expected to occur by chance [63]. It ranges from −1 to +1, where +1 indicates perfect agreement and values of zero or less indicate a performance no better than random [63]. While the k is recognized for its simplicity, the fact that both commission and omission errors are accounted for in one parameter, and its relative tolerance to zero values in the confusion matrix, it is also criticized for being inherently dependent on prevalence and that this dependency introduces bias and statistical artefacts to estimates of accuracy [62].

Ensemble Model Development
The TSS has been shown to provide a more robust measure of performance than the other three metrics that we calculated, particularly the AUC [25,59]. Given that the method to select algorithms to create an ensemble is still a gray area, we chose to use a metric that has been shown to work better in previous studies [22]. We constructed ensemble maps from all the 17 algorithms for each month, using the TSS as a weighting factor [22]. This approach gave higher statistical weights to algorithms whose TSS values were high, and lower weights to algorithms whose TSS values were low. The analysis workflow is illustrated in Figure 2. specificity, and is insensitive to prevalence [62]. Like kappa, it takes into account both omission and commission errors, and success as a result of random guessing, and ranges from −1 to +1, where +1 indicates perfect agreement and values of zero or less indicate a performance no better than random. [62]. The Cohen's kappa is a widely used statistic which corrects the overall accuracy of model predictions by the accuracy expected to occur by chance [63]. It ranges from −1 to +1, where +1 indicates perfect agreement and values of zero or less indicate a performance no better than random [63]. While the k is recognized for its simplicity, the fact that both commission and omission errors are accounted for in one parameter, and its relative tolerance to zero values in the confusion matrix, it is also criticized for being inherently dependent on prevalence and that this dependency introduces bias and statistical artefacts to estimates of accuracy [62].

Ensemble Model Development
The TSS has been shown to provide a more robust measure of performance than the other three metrics that we calculated, particularly the AUC [25], [59]. Given that the method to select algorithms to create an ensemble is still a gray area, we chose to use a metric that has been shown to work better in previous studies [22]. We constructed ensemble maps from all the 17 algorithms for each month, using the TSS as a weighting factor [22]. This approach gave higher statistical weights to algorithms whose TSS values were high, and lower weights to algorithms whose TSS values were low. The analysis workflow is illustrated in Figure 2.

Results
Our averaged model performance results for the monthly models for March to November, constructed using 17 algorithms are shown in Figure 3

Results
Our averaged model performance results for the monthly models for March to November, constructed using 17 algorithms are shown in Figure 3 Table S1, and the TSS values are in bold print. The mean ROC plots for the 17 algorithms in each of the 9 monthly models are shown in Figure S1, with confidence limits within 2 standard deviations of the mean. The ROC plots are probability plots of sensitivity (true positive rate) against 1-specificity (false positive rate) where an algorithm which presents a curve closer to the top-left corner indicates better performance compared to one where the curve is closer to the 45-degree line of the ROC space. For instance, the July ROC plots generated with the generalized linear model (GLM), FDA and MAXL ( Figure S1:(5,68,104)) algorithms indicate ROC plots whose part of the plot lean towards the 45-degree line, hence lower performance. These results are quite different from the MAXENT ROC plots ( Figure S1:(145-153)) for all the months, which illustrate higher performance. The variable response curves for SSC, SSHA and SST are shown in Figures S2, S3 and S4 respectively, where the values which influenced model performances are indicated by the peaks of each curve in a specific algorithm. However, there are algorithms which show relatively monotonic shapes for certain variables, e.g., the BRT and MAXL for the SSC, SSHA and SST variables. The shape, peak, and range of values around the peak present important information which influences the HSI map. Figure S5 illustrates the contribution of each variable in the respective model. These results illustrate the relative importance of each of the 3 variables used to construct the model, hence variables with high values on a plot indicate a relatively higher contribution by that variable to the model fit and HSI map, and vice-versa. Notable, for instance, is the importance of SSHA in the GLM, SVM, MLP and RPART algorithms between April and June, while SST shows a strong variable contribution in MAHA, MAXENT, MAXL, and BIOC between September and November. The ensemble HSI maps for March to November period are shown in Figure 4, which illustrates the aggregated map. The high value pixels (HSI>0.7) show a consistent pattern where the latitudinal displacement increases northward from March toward October and November.

Discussion
We modelled skipjack tuna habitats in the western North Pacific from March to November, using fishery data and satellite remotely sensed variables in 17 algorithms implemented in the sdm package [15]. Our approach employed ensemble modelling techniques drawn from correlative and machine-learning algorithms to generate an ensemble habitat prediction. This process enabled us to explore the inherent characteristics of each algorithm as well as the ensemble modeling approach for applications such as potential fishing zone forecasting and fishery management options. Given the collection of algorithms we applied, our work explores a much wider scope of habitat algorithms for skipjack tuna in this area than previously. Ensemble modelling has evolved from the need to aggregate outputs from different models, optimize the success rates while minimizing the variations and inherent modelling limitations of constituent models. Since models are representations of reality, and no single model can be perfect in estimating real world conditions, using a number of models is deemed an appropriate way of optimizing the true "signal" about the relationships a model is aiming to capture, and minimizing some "noise" created by errors and uncertainties in the data and the model structure [64]. Generally, consensus on the ideal method to aggregate outputs of single algorithms into an ensemble is still lacking. Some authors find simple averages adequate [64] while others use a threshold to select better performing algorithms first, followed by weighting of these map outputs into an ensemble [22], [64], [65]. Indeed, the intended use of an ensemble prediction may largely determine the aggregating method. Nevertheless, the practical applications (e.g., fishery and marine conservation) of ensemble model outputs are clear, and include reducing the likelihood of making decisions based on maps that are far from the 'truth', the costs of relying on a single

Discussion
We modelled skipjack tuna habitats in the western North Pacific from March to November, using fishery data and satellite remotely sensed variables in 17 algorithms implemented in the sdm package [15]. Our approach employed ensemble modelling techniques drawn from correlative and machine-learning algorithms to generate an ensemble habitat prediction. This process enabled us to explore the inherent characteristics of each algorithm as well as the ensemble modeling approach for applications such as potential fishing zone forecasting and fishery management options. Given the collection of algorithms we applied, our work explores a much wider scope of habitat algorithms for skipjack tuna in this area than previously. Ensemble modelling has evolved from the need to aggregate outputs from different models, optimize the success rates while minimizing the variations and inherent modelling limitations of constituent models. Since models are representations of reality, and no single model can be perfect in estimating real world conditions, using a number of models is deemed an appropriate way of optimizing the true "signal" about the relationships a model is aiming to capture, and minimizing some "noise" created by errors and uncertainties in the data and the model structure [64]. Generally, consensus on the ideal method to aggregate outputs of single algorithms into an ensemble is still lacking. Some authors find simple averages adequate [64] while others use a threshold to select better performing algorithms first, followed by weighting of these map outputs into an ensemble [22,64,65]. Indeed, the intended use of an ensemble prediction may largely determine the aggregating method. Nevertheless, the practical applications (e.g., fishery and marine conservation) of ensemble model outputs are clear, and include reducing the likelihood of making decisions based on maps that are far from the 'truth', the costs of relying on a single algorithm prediction compared to an ensemble, and the benefits of using ensembles to explore the impacts of future phenomena such as fishery forecasts and impacts of climate change [66].
Our results showed that the performance of the 17 single algorithm models varied, with SVM, BRT, RF, MARS, GAM, CART, MLP, RPART, and MAXENT, showing consistently high performance in all the months (Figure 3). However, we also noted lower performance for FDA, MDA, BIOC, DOM, MAXL, MAHA and RBF algorithms in a number of months, particularly in May, July and August, October and November (Figure 3). Algorithm performance has been the subject of some studies, where some machine learning algorithms such as MAXENT and RF have been shown to outperform others [61,67]. In this work, four performance metrics were employed to measure algorithm performance, with the understanding that each metric measures different aspects of performance [67]. The COR metric for the MAHA algorithm was low in all months. However, previous work has shown that the TSS is a more reliable metric than the AUC, COR, and Kappa, because it is independent of prevalence and reflects true ecological phenomena rather than statistical artefacts [59,62], hence we relied on the TSS to generate ensemble predictions. Algorithm performance can be influenced by a number of factors, namely the inherent mathematical configuration, the environmental distribution of the data in the region of interest, in the training sample, and in regions that might be used for projection, as well as the input parameter settings [67]. For instance, it has been observed that MAXENT's better performance in comparison to the other algorithms might be partly due to how the environmental variables and their interactions are modelled, where more mathematical complexity of the model is progressively improved when more data are available [68,69]. In other studies, it has been pointed out that generative methods such as MAXENT and RF render better results with small sample sizes, maybe due to faster convergence to their higher asymptotic error than discriminative methods [68]. By contrast, discriminative methods such as GLM and GAM improve their accuracy as the number of records increases and may even surpass results offered by generative methods at large sample sizes [68,70]. In this study, we provided all algorithms with the same input datasets and configurations as implemented in the sdm package and thus assessed the differences in output, which we largely attribute to the individual algorithm behavior. It is also worth noting that our results depict outputs of the currently implemented versions of the 17 algorithms as implemented in R, and it is likely that future updates of the algorithms may lead to variations in observed performance.
The response curves and variable importance plots ( Figures S2, S3, S4 and S5) showed consistencies by some algorithms on the importance of either SST, SSC or SSHA over time. For instance, the GLM, SVM, MLP and RPART models showed that SSHA had a higher importance between March and June, while the RBF, MAHA, MAXL, DOM, BIOC and FDA placed higher variable importance on SST in May and June. By contrast, between September and November, a number of algorithms placed higher variable importance on SST and SSC. Visualizing the fitted functions of the various variables was important in exploring the modelled relationships, and also understanding the differences among the methods and how different fitted functions influenced mapped predictions [67]. From an ecological perspective, between March and June, skipjack tuna are abundant in the Kuroshio area, which is dominated by warm waters, as they advance north during the northward migration, hence the models may not have been sensitive to the SST signal as an important variable during this period, an observation made by earlier work [5]. However, as the fish migrated into the Kuroshio Oyashio Transition Zone (June-August), and later towards the Oyashio area (September-November) where the waters are cooler, the importance of SST was more apparent [8]. Skipjack tuna are warm water species, and the mechanisms with which they forage within warm waters along productive thermal and ocean color fronts are well explained in previous work [29,30,36,71]. Furthermore, it is worth noting that some algorithms did not do well at detecting responses on some variables, hence displayed monotonic response functions. For instance, the BRT and CART (March to June) algorithms showed a flat response for SSC in all months, yet other algorithms were able to pick that signal. A similar response was noted with the BRT algorithm for April, June, October and November, as well as the CART and MAXL algorithms between April and June. In another study which compared the performance of 6 algorithms (5 of which we also used), it was noted that individual models provided slightly different results with the same dataset [25]. This underscores the importance of multi-algorithm approaches, compared to relying entirely on single algorithm models [9,25].
The ensemble predictions ( Figure 4) were created from a suite of 17 species distribution algorithms which were weighted with the TSS values (Table S1). We used the TSS weighting for creating ensemble predictions since there isn't a clearly defined selection criteria and weighting metric for creating ensembles [22]. The TSS metric has been shown to be more reliable than the AUC, and to provide a more reasonable output in other studies [22,25]. In addition to these considerations, we also concur with [67] that the application for which an ensemble model is sought constitutes a fundamental question that shapes the choice and weighting of the single algorithms that create the ensemble prediction. The ensemble prediction maps aggregate the outputs of the individual algorithms, and are comparable in pattern to maps derived from previous work in the same area [5]. Furthermore, the observed northward displacement of high suitability indices from March to October is also consistent with recent descriptions of skipjack tuna migrations in the same area [36]. However, comparisons of ensemble outputs in a multi-year study could add more information on interannual variability of habitats. The foreseeable applications of our findings in the western North Pacific include potential fishing zone modelling and forecasting, and prediction of range shifts driven by climate change. Potential fishing zone forecasting comprises a much shorter time scale compared to climate change driven range shifts, but for both, it is important to understand how the algorithms perform when projected into new environmental spaces not sampled by the training data [67].
Our work was subject to a number of data and modelling limitations which we point out here for future consideration. First, we used a single-year fishery dataset as a response variable due to logistics related to data availability for this study. An inter-annual dataset using a multi-algorithm approach would have been more useful in elucidating the strengths and limitations of each algorithm, as well as the inter-annual variability in habitat characteristics as depicted by each algorithm. This would also have provided an opportunity for validations using an independent dataset as recommended by [9]. In addition, like any other field dataset, fishery-derived data are subject to sampling biases associated with sampling effort distribution, either influenced by fisherfolk choices or weather considerations, or both [26]. Second, our study used a set of 3 satellite-derived predictor variables (SST, SSC, SSHA) which we acknowledge are not the only environmental variables that are important for skipjack tuna habitat characterization. Indeed, other studies have also used variables such as oxygen concentration, currents, sub-surface temperatures and mixed layer depths [22,23,36]. We limited our work to only SSC, SST, SSHA since these are widely applied predictors in tuna habitat studies [7,72], which also enabled us manage the algorithm and predictor variable permutations in our runs. Third, all the response and predictor variable layers were standardized to~25-km spatial resolution, monthly averaged files to create consistency in spatial-temporal resolutions, while cognizant of the fact that resolution trade-offs are often inevitable in species distribution modelling [26]. Therefore, it is possible to use a higher resolution dataset, especially for smaller areas where there is need to resolve fine-scale oceanographic features. Fourth, there are biases introduced by mismatches in space and time between the sampling protocols for the response and predictor variable datasets [20]. For example, while the fishery data were recorded at point spatial resolution, on a daily basis, for the purposes of this study, they were aggregated to~25-km resolution to match the resolutions of the environmental predictors. At times, while temporal averaging to monthly timescales is deemed necessary to improve data coverage due to cloud cover (e.g., for SSC), it may hinder a model's ability to resolve fine-scale daily to weekly signals that influence habitat characterization [5,26].

Conclusions
This work is a first attempt at a multi-algorithm species distribution habitat modelling for skipjack tuna in the western North Pacific, using presence-only response data and satellite-derived predictor variables. We draw the following conclusions within the confines of the datasets we used, the algorithms applied, and the input parameter settings of the algorithms applied in the study area. First, the multi-algorithm approach elucidated the strengths and weaknesses of the various algorithms to model and predict skipjack tuna habitats. According to the performance statistics, the SVM, BRT, RF, MARS, GAM, CART, MLP, RPART, and MAXENT achieved consistently better performance than other algorithms. Some of the lowest-performing algorithms were the FDA,