Spatiotemporal Patterns in the Distribution of Albacore, Bigeye, Skipjack, and Yellowﬁn Tuna Species within the Exclusive Economic Zones of Tonga for the Years 2002 to 2018

: The Tongan ﬁsheries targeting the species of albacore


Introduction
Tonga is enveloped within the geographical boundary of the Western Central Pacific Ocean (WCPO) and its exclusive economic zone (EEZ) is managed by its Ministry of Fisheries with the objective of maintaining the viability of each tuna fishery through the implementation of its management policies.The most economically important tuna species in Tonga are albacore, bigeye, skipjack, and yellowfin, which account for over 95% of all tuna fisheries by annual catch and economic value.Tuna production is the largest commercial fishery in Tonga and is estimated at 2000 metric tonnes (mt) per year (approximately 17% of 30% of Tonga's marine resource-related benefits) [1].Longline is the main method for commercial tuna fishing, starting from the 1970s, although pole and line has also been used [1,2].
The application of satellite remote sensing oceanographic data within the geographic information system (GIS) has been identified as a tool that may improve tuna production and management outcomes [3,4].Potential fishing zones (PFZ) have been identified in the EEZ of Madagascar, the Arabian Sea, Tropical Pacific Ocean, the north-west coast of India, and the western North Pacific Ocean, using satellite data collected from environmental parameters and oceanographic conditions [3,5,6].As stated above, biophysical environmental parameters such as sea surface temperature (SST), chlorophyll-a (Chl-a) concentration, and water depths play an important role in influencing the spatial and temporal distributions of tuna [4,7].Predictions of PFZs are established based on species distribution models aimed at parameterizing the environmental factors that influence the biology and habitats of the species of interest.
Studies have shown that environmental variables and oceanic conditions influence tuna presence; thus, characterizing their habitats is essential to understanding their spatial and temporal distributions [8][9][10].Sea surface temperature and Chl-a concentration have been identified as playing a role in generating spatial and temporal distribution patterns of tuna in the WCPO [11], East Pacific Ocean [12], and the East Coast of Australia [7].These studies have indicated that tropical and subtropical tuna prefer warm waters and are found in regions with high Chl-a concentration.In addition, it has been shown that the seasonal and annual patterns of tuna abundance are affected by SST and Chl-a and that their variability is related to ocean upwelling and El Nino events [13][14][15].Utilization of satellite remotely sensed and fisheries data for management and assessment of tuna and using the information as basis for decision support by managers has not been carried out in Tonga [1] nor in its surrounding island countries like Fiji, Samoa, and Niue.
Chlorophyll-a plays a crucial role in determining the distribution of tuna as it serves as a primary indicator of primary productivity in the ocean, influencing the abundance and distribution of their prey species [13,16].Hence, tuna rely on the presence of chlorophyll-arich areas as they indicate the presence of productive feeding grounds, which ultimately affects their migration patterns and distribution in search of optimal foraging opportunities.Sea surface chlorophyll-a gradients (chl_grad), affected by nutrient and light availability as well as oceanic currents, offer valuable insights into the distribution, abundance, and health of marine primary producers, contributing to the assessment of ecosystem productivity [17].Sea surface temperature significantly influences tuna distribution as these species are ectothermic, meaning their body temperature is dependent on the surrounding environment [18].Optimal temperature ranges dictate the presence and movement of tuna, with warmer waters often associated with increased abundance, as they provide favorable conditions for their metabolic activity and prey availability [19][20][21].Conversely, deviations from the preferred temperature range can lead to shifts in tuna distribution, impacting their migration patterns and potentially reducing their availability in certain regions [13,20].Ocean temperature in Tonga is above 18 • C and is favorable for tropical and subtropical tuna species, which are known to have temperature preference range of ≥25 • C [20].Albacore, which is known to be temperate species, is also caught in Tongan waters and is known to have temperature preference range between 20 • C and 30 • C in each hemisphere [20].The sea surface temperature gradient (sst_grad) quantifies the speed of temperature change in the ocean over a specific distance, with higher gradients indicating faster changes and lower gradients indicating more gradual changes [22].These gradients are influenced by factors such as solar radiation, ocean currents, and atmospheric conditions, and they significantly impact weather patterns, ocean circulation, and the distribution of marine life [23].Bathymetry has been identified to affect catch rate of yellowfin and bigeye tuna in the Indian Ocean, with their catch per unit effort (CPUE) generally becoming higher as steep bathymetry zones become larger [24].Hence, bathymetry may also play a role as Tonga EEZ envelops the Tonga Trench, Tonga Ridge, Tofua Arce Volcanic Front, the northern end of the Tonga Kermadec Arc and the westward region of the Lau Basin.These geologically bathypelagic features are part of the island nation's fishing grounds and may influence its oceanic conditions such as; surface water temperature, nutrient, salt content, and mixed layer depth, which are habitats for pelagic species and fishing fleets.Bathymetry may also play a role since volcanic seamount lines run through the fishing ground and include very deep water more than 5000 m deep and very shallow seamounts and a large shelf that is about 2000 m deep.
Apex predators in the marine environment, including tuna, are shown to exploit frontal zones; thus, for both SST and Chl-a, frontal gradients were computed and used to model the distribution of the four tuna species [25,26].Although the range of species distribution models can be used for the purpose of the current project, generalized additive models (GAMs) were used instead.GAMs allow for non-linear relationship between response variables, in this case CPUE of population of tuna species, and sets of predictors without explicitly specifying the form of the relationship, thus enhance our knowledge of the influence of environmental variables on the distribution and abundance of populations [27,28].They are semiparametric extensions of generalized liner models with an assumption that the functions are additive and the components are smooth [29].GAMs have the ability to handle highly nonlinear and non-monotomic relationships between responses and explanatory variables, making them ideal for expressing underlying relationships in ecological systems [19].This technique has been successfully used to describe the habitat preferences of target and non-target species and has often been used to predict species distributions [30].The main objectives of this study are to identify (1) the spatial and inter-annual variability in the CPUE of the four species within the EEZ of Tonga and (2) the form and strength of influence of the above mentioned environmental variables.

Study Region
This study was conducted in the archipelago waters of Tonga, a small island country located 700 km south-west of Fiji and 1900 km north-west of New Zealand, in the South Pacific Ocean.It covers an EEZ (latitude 14.5 • S-20.22 • S, longitude 171.31 • W-179.10 • W) of approximately 800 km 2 (Figure 1).This EEZ envelops the northern end of the Tonga Trench, Tonga Ridge, Tofua Arce Volcanic Front, the northern end of the Tonga Kermadec Arc, and the westward region of the Lau Basin.Two geologically parallel north-to-south chains of volcanic seamounts run along the Tonga Ridge and the famous seamount of Capricorn about 193 km east of Vava'u Island.This fishing ground supports the nation's commercial tuna fishery harvest of about two thousand tons per year.

Fishery Data
Datasets for albacore, bigeye, skipjack, and yellowfin tuna catch from January 2002 to December 2018 were used to model the distribution of the four species in the EEZ of Tonga.The catch and effort data were compiled by the Tongan longline fishery and provided by the Tonga Ministry of Fishery and the South Pacific Community (SPC) Office in New Caledonia.The entire fish catch was checked by a minimum of two fisheries offices to ensure compliance [1].The fishery data of the 1 • spatial grid comprising daily fishing positions (latitude and longitude), fishing effort (in number of hooks), fishing date (in day, month, and year), and catch (in numbers).The observed catch rate was calculated as the weight of the catch in metric tons divided by the number of hooks deployed per fishing record, providing a standardized measure of fishing efficiency and effort in capturing the target species [31,32].The CPUE data were then aggregated into the monthly level to match the temporal scales of the environmental data.The distributions of the CPUEs of albacore, bigeye, skipjack, and yellowfin are shown in Figure 2, respectively.

Fishery Data
Datasets for albacore, bigeye, skipjack, and yellowfin tuna catch from January 2002 to December 2018 were used to model the distribution of the four species in the EEZ of Tonga.The catch and effort data were compiled by the Tongan longline fishery and provided by the Tonga Ministry of Fishery and the South Pacific Community (SPC) Office in New Caledonia.The entire fish catch was checked by a minimum of two fisheries offices to ensure compliance [1].The fishery data of the 1° spatial grid comprising daily fishing positions (latitude and longitude), fishing effort (in number of hooks), fishing date (in day, month, and year), and catch (in numbers).The observed catch rate was calculated as the weight of the catch in metric tons divided by the number of hooks deployed per fishing record, providing a standardized measure of fishing efficiency and effort in capturing the target species [31,32].The CPUE data were then aggregated into the monthly level to

Environmental Data
To assess the influence of environmental variables on the distribution of the four tuna species, SST and Chl-a were downloaded from the NOAA repository [33].Both of these data sources are at a monthly resolution.In addition, to assess the potential influence of fronts regions, frontal gradients for SST and Chl-a were computed for each monthly layer, using the grec R package, version 1.5.0 [34].These were then normalized by dividing the values by the maximum for the month.The distributions of these environmental variables are shown in Figure 3, respectively.These data made it possible to identify well-defined fronts [35].Before applying distribution models, the existence and strength of multi-collinearity were assessed.This was completed by computing the variance inflation factor (VIF).There was no indication of the existence of serious multi-collinearity as none of the computed VIF values were more than 5 (V IF ≤ 5).Seven GAM models were considered to model the distribution of the four species.Selection among the seven models was based on k-fold cross-validation (where k is repeated 5 times), using root mean squared prediction error (RMSPE) as a performance measure.The seven models considered were:

Environmental Data
To assess the influence of environmental variables on the distribution of the four tuna species, SST and Chl-a were downloaded from the NOAA repository [33].Both of these data sources are at a monthly resolution.In addition, to assess the potential influence of fronts regions, frontal gradients for SST and Chl-a were computed for each monthly layer, using the grec R package, version 1.5.0 [34].These were then normalized by dividing the values by the maximum for the month.The distributions of these environmental variables are shown in Figure 3, respectively.These data made it possible to identify well-defined fronts [35].Before applying distribution models, the existence and strength of multi-collinearity were assessed.This was completed by computing the variance inflation factor (VIF).There was no indication of the existence of serious multi-collinearity as none of the computed VIF values were more than 5 ( ≤ 5).Seven GAM models were considered to model the distribution of the four species.Selection among the seven models was based

Index of Variable Importance
There are multiple approaches used to compute the relative importance of variables.One such approach applied in the context of species distribution modelling is implemented in the R packages SSDM [36], biomod2 [37].The approach determines the relative importance of each variable by comparing model predictions under the full model (model with all the variables) to that based on model without the variable of interest.The change in correlation between predictions under the full model to that after removing a predictor indicates how important the removed variable was.The index of importance was computed as shown below: where    is the relative index of importance,   prediction from the best model,    is the prediction from the variable i permuted.To obtain a measure of variability of the relative importance required repeating this process 50 times.All the processing, analysis, visualization, and report generation were performed using multiple packages [38][39][40][41][42][43][44][45][46] in R software [47].

Index of Variable Importance
There are multiple approaches used to compute the relative importance of variables.One such approach applied in the context of species distribution modelling is implemented in the R packages SSDM [36], biomod2 [37].The approach determines the relative importance of each variable by comparing model predictions under the full model (model with all the variables) to that based on model without the variable of interest.The change in correlation between predictions under the full model to that after removing a predictor indicates how important the removed variable was.The index of importance was computed as shown below: where I r i is the relative index of importance, P r prediction from the best model, P p i is the prediction from the variable ii permuted.To obtain a measure of variability of the relative importance required repeating this process 50 times.All the processing, analysis, visualization, and report generation were performed using multiple packages [38][39][40][41][42][43][44][45][46] in R software [47].

Exploration of the Occurrence
An overview of the distribution of CPUE used in all the distribution modelling is provided in Figures S1-S4 for albacore, bigeye, skipjack, and yellowfin, respectively.The first available data for NASA Aqua was in June 2002, hence the overlay of catch data started over the Tonga EEZ.The figures show monthly and annual distributions of the CPUE of the four species by the longline fishery throughout the year from January to December.The distribution of CPUE shown suggests that most of the fishing effort of the tuna longliners was concentrated in the 15-24 • S. Bigeye, skipjack, and yellowfin CPUE shows a central-northernmost distribution, and was caught mainly between 14 • S and 22 • S. The albacore species thrives in temperate to subtropical regions and exhibits its highest catch rates between the central-southern part of the EEZ.The highest CPUE for all species are observed from 2002-2008 and from 2013-2018 and the lowest CPUE are observed from 2009-2012.In addition, the figure also shows a broad scale seasonality in the CPUE of the four tuna species and an increase in their CPUE especially over the recent period.

Exploration of Environmental Variables
The spatiotemporal pattern in SST, Chl-a, sea surface temperature gradient, sea surface chlorophyll gradient, and bathymetry data within the EEZ of Tonga are shown in Figures S5-S9, respectively.The distribution of SST is not homogeneous throughout the EEZ of Tonga, with pronounced seasonal variation.The value of the temperature ranges from 18 • C to 30 • C. From December to April, the temperature remains between 27 • C and 30 • C. From May, it falls significantly to reach its minimum in August, at approximately 18 • C. SST rises from September until November.The northern part of the EEZ is warmer compared to the central and southern parts.The local relative abundance of the four tuna species increases towards the central part of the EEZ (Figures S1-S4), indicating a preference of tropical tuna species to warmer waters.The distribution of SST shows that, during cooler months (April-October), the fishing operations aggregated in more temperate waters of SST 22 • C to 30 • C, where the highest catches occurred.The value of Chl-a concentration was at a range of 0.05-0.1 mgm −3 and was observed to be homogeneous and with no clear spatial and temporal patterns.However, generally there was a higher Chl-a concentration in the central part of the EEZ during summer, from November to January.
Bathymetry in this fishing ground is in the range of approximately 2000−8000 m.The islands of Tonga geographically distributed along the south-west to north-east of the EEZ.Very deep water of more than 5000 m covers the Tonga Trench and to the south-eastern part of the EEZ, with very shallow volcanic sea mounts and a large shelf ≥2000 m deep at northwestern part.The sea surface chlorophyll-a gradient ranged from 0 to 0.75 mgm −3 m −1 , with slightly higher gradients occurring in the summer months (November to March) and an increase observed during the cooler months (May to August) in the most recent period of the study.The sea surface temperature gradient exhibited a wide annual variation ranging from 0 to 0.60 • Cm −1 with slightly higher values observed in the periods of 2004-2010 and 2014-2017.

Model Performance
The performance of each of the seven models, from k-fold cross-validation as measured by the RMSPE, is presented in Table 1.Most models with only environmental variables and bathymetry as predictors appear to have comparable performance with noticeable improvement in performance with the addition of catch month and year (e.g., from 11.703 (model 5) to 10.293 (model 7) for albacore).Figure 4 shows the partial effects of predictors on the CPUE of all the four tuna stocks.The partial effects with their partial residuals are presented in Figure S10.For ease of flow, partial effects of all the predictors are presented together.The partial effects of year indicated a gradual increase in CPUE from 2002 to 2018.CPUE appeared to decrease with increases in SST from 20 • C to about 25 • C, after which it appeared to increase again.Relatively higher CPUE were also associated with low SST fronts.Albacore appeared to be associated with relatively low chlorophyll water and also a relatively stable chlorophyll front.Within the EEZ of Tonga, albacore appeared to be associated with areas that are deeper with CPUE declining in shallower water <2000 m.Seasonally, higher CPUE were observed around the middle of calendar year.For bigeye tuna, the partial effects of year indicates a gradual increase in CPUE from 2002 to 2007, after which it declined and only started to increase after 2015.The CPUE of bigeye tuna appears to increase with increases in SST up to 25 • C, after which it remained stable.Relatively higher CPUE were also associated with low SST fronts.Bigeye tuna appear to be associated with relatively high-chlorophyll water and also a relatively stable chlorophyll front.Within the EEZ of Tonga, bigeye appear to be associated with areas that are deeper, with CPUE declining linearly toward shallower water.Seasonally higher CPUE were observed around the start of the calendar year.Plots of predicted values vs. observed values for each species presented in log-scale are shown in Figure S11.
Table 1.Performance of the seven models considered for albacore, bigeye, skipjack, and yellowfin in the EEZ of Tonga.The mean RMSE and mae (mean absolute error) of the test data are presented.The partial effects of predictors for skipjack tuna are presented in Figure 4.The partial effects of year indicate a gradual decrease in CPUE from 2002 to 2008, after which it gradually increased.The CPUE of skipjack tuna appears to decrease with increases in SST up to 25 • C, after which it remained relatively stable.The partial effects of chlorophyll indicates that the CPUE of skipjack gradually declined with an increase in log Chl.Although the effects of chl_grad are largely flat, they suggests an association of higher CPUE with a low and strong chlorophyll frontal area.A similar pattern appears to characterize the partial effect of sst_grad on the CPUE of skipjack.Within the EEZ of Tonga, skipjack tuna appear to be associated with areas that are around 5000 m deep with CPUE declining on both sides.Seasonally higher CPUE were observed around July and August.C, after which it declined.However, sst frontal (sst_grad) appears to have limited/minor influence on the CPUE of yellowfin tuna.A higher CPUE of yellowfin tuna was associated with water mass characterized by strong chlorophyll fronts (chl_grad).In addition, the CPUE of yellowfin tuna also appears to increase with increases in chlorophyll concentration (logChl).Within the EEZ of Tonga, the CPUE of yellowfin tuna appears to remain stable up to a depth of 5000 m, declining in waters deeper than 5000 m.Seasonally, higher CPUE were observed around July and August.

Species
The index of relative importance of the predictors influencing the distribution and abundance of the four tuna populations caught within the EEZ of Tonga is shown in Figure 5.The numerical index of relative importance of predictors for the four tuna species is presented in Table S1.As can be seen in the figure, catch year and catch month were the two most important predictors for all the three species, except for albacore tuna, for which the order reverses.The third most important variable was bottom depth for skipjack and yellowfin tuna and SST for albacore and bigeye tuna.The remaining variables appear to have a minor role.Model diagnostics for the best model were also checked and there was no substantial violation of the assumption of independent and identically distributed variables.Plots of linear predictor vs. deviance residuals and quantile-quantile of residuals are shown in Figure S12 for all species.Figure 6 shows plots of a posterior predictive distribution of the models for the four tuna species.It illustrates that, despite the relatively low deviance explained (Tables 2 and S2) for all the tuna populations modelled, the model captures the dynamic of the data generation process.Table 2 presents the cumulative proportion of the deviance explained from the null model to the full model.Table S2 presents smooth effects and Table S3 presents the parametric effects of the best GAM models fitted for each tuna species.As is evident from Table 2, the full models (models including all terms) demonstrate the highest cumulative deviance explained across all species, with albacore having the highest value at 27.622, followed by skipjack and yellowfin, while bigeye has the lowest value at 12.183.Also, it can be seen (Table S2) that all of the smooth terms were significant for albacore (with p-value < 0.05).For bigeye tuna, sst_grad was not statistically significant, whereas the remaining terms had a significant effect.The effects of both frontal terms (sst_grad and chl_grad) and chlorophyll were not significant for skipjack tuna with SST, bottom depth, and month having statistically significant effects.For yellowfin tuna, all the terms, except for sst_grad, had a significant effect on its spatiotemporal pattern.

Standardized Index of Abundance
The standardized index of abundance, based on the best GAM model, is presented in Figure 7.For comparison, both the standardized and nominal index of abundance were normalized, I n = I/ I .The trends for both albacore and bigeye tuna are characterized by inter-annual variability with no long-term trend over the period considered.Skipjack and yellowfin tuna appear to show an overall increase in abundance over the period 2002-2018.

Map of Predicted Suitable Habitat
Maps of predicted suitable habitat for albacore, bigeye, skipjack, and yellowfin are shown in Figures 8-11.They are arranged quarterly to illustrate how their abundance changes with season.Albacore tuna appeared to be characterized by seasonal variability in CPUE where they become more abundant during the months of July to August, with regions in the northern and southern parts of the EEZ characterized by relatively higher abundance (Figure 8).Conversely, as shown in Figure 9, bigeye tuna appeared to be more abundant during the first half of the year (January-June) and in regions in the southern and eastern part of the EEZ.Although skipjack increased in abundance after 2013 (Figure 10), they appeared to be characterized by the same seasonal cycle as albacore tuna and peak in abundance around July to August, but are associated with the same regions as bigeye tuna.Yellowfin tuna has also been increasing in abundance since about 2010; although characterized by inter-annual variability, they appear to be in higher abundance from December to July and mostly around the island chain (Figure 11).

Map of Predicted Suitable Habitat
Maps of predicted suitable habitat for albacore, bigeye, skipjack, and yellowfin are shown in Figures 8-11.They are arranged quarterly to illustrate how their abundance changes with season.Albacore tuna appeared to be characterized by seasonal variability in CPUE where they become more abundant during the months of July to August, with regions in the northern and southern parts of the EEZ characterized by relatively higher abundance (Figure 8).Conversely, as shown in Figure 9, bigeye tuna appeared to be more abundant during the first half of the year (January-June) and in regions in the southern and eastern part of the EEZ.Although skipjack increased in abundance after 2013 (Figure 10), they appeared to be characterized by the same seasonal cycle as albacore tuna and peak in abundance around July to August, but are associated with the same regions as bigeye tuna.Yellowfin tuna has also been increasing in abundance since about 2010; although characterized by inter-annual variability, they appear to be in higher abundance from December to July and mostly around the island chain (Figure 11).

Discussion
Our purpose in modeling albacore, bigeye, yellowfin, and skipjack distribution in the EEZ of Tonga was to link tuna CPUE to environmental, physical, and temporal factors and subsequently to generate predictive habitat suitability maps and standardized indices of abundance, which could be useful for national planning and policy development for conservation and sustainable management of the species.The environmental variables included in this study were SST, Chl-a, sst grad, chl_grad and bottom depth (as physical predictors).Other variables (e.g., salinity levels, mixed layer depth, sea surface height, ocean currents, dissolved oxygen, and others) have been employed in previous studies as predictors when modeling the suitability of tuna habitats [49][50][51][52][53][54].Studies indicate that these surface and subsurface layers of oceanic environments also influence the distribution and presence of tuna in the water column.The movement of ocean surface waters, including converging and diverging flows, is linked to fluctuations in sea surface height.Tuna utilize these variations as signals to identify areas where food is plentiful.[49].The movement of ocean currents has an effect on the mobility of phytoplankton and small pelagic fish, thus shaping the dispersion of marine organisms at higher trophic levels [50].The distribution and presence of tuna in the ocean environment are also influenced by dissolved oxygen levels [51] and salinity [52,53].However, when accessible, these investigations are constrained either to a single species or to relatively small regions [7].Furthermore, when these variables are integrated with sea surface temperature and chlorophyll concentration to model habitat suitability, research has indicated that SST and Chl-a concentration are the predominant factors among them [3,10,19,[55][56][57].Our study also included extensive tuna longline datasets with spatiotemporal coverages, where the information was used to model tuna habitat preferences and distribution as a number of studies have shown [58][59][60].It may be premature to draw strong inferences regarding environmental effects on tuna distribution and abundance, even so, our results indicated that temporal variables, e.g., inter-annual variability and seasonality are influential.Bottom depth served as surrogate for environmental conditions not captured by SST, Chl-a, and the two frontal indices (sst_grad, and chl_grad).This played an important role following year and month for all the tuna species considered except for albacore tuna.SST was the most important variable following year, month, and bottom depth for all the tuna species considered in this study.

Discussion
Our purpose in modeling albacore, bigeye, yellowfin, and skipjack distribution in the EEZ of Tonga was to link tuna CPUE to environmental, physical, and temporal factors and subsequently to generate predictive habitat suitability maps and standardized indices of abundance, which could be useful for national planning and policy development for conservation and sustainable management of the species.The environmental variables included in this study were SST, Chl-a, sst grad, chl_grad and bottom depth (as physical predictors).Other variables (e.g., salinity levels, mixed layer depth, sea surface height, ocean currents, dissolved oxygen, and others) have been employed in previous studies as predictors when modeling the suitability of tuna habitats [49][50][51][52][53][54].Studies indicate that these surface and subsurface layers of oceanic environments also influence the distribution and presence of tuna in the water column.The movement of ocean surface waters, including converging and diverging flows, is linked to fluctuations in sea surface height.Tuna utilize these variations as signals to identify areas where food is plentiful.[49].The movement of ocean currents has an effect on the mobility of phytoplankton and small pelagic fish, thus shaping the dispersion of marine organisms at higher trophic levels [50].The distribution and presence of tuna in the ocean environment are also influenced by dissolved oxygen levels [51] and salinity [52,53].However, when accessible, these investigations are constrained either to a single species or to relatively small regions [7].Furthermore, when these variables are integrated with sea surface temperature and chlorophyll concentration to model habitat suitability, research has indicated that SST and Chl-a concentration are the predominant factors among them [3,10,19,[55][56][57].Our study also included extensive tuna longline datasets with spatiotemporal coverages, where the information was used to model tuna habitat preferences and distribution as a number of studies have shown [58][59][60].It may be premature to draw strong inferences regarding environmental effects on tuna distribution and abundance, even so, our results indicated that temporal variables, e.g., inter-annual variability and seasonality are influential.Bottom depth served as surrogate for environmental conditions not captured by SST, Chl-a, and the two frontal indices (sst_grad, and chl_grad).This played an important role following year and month for all the tuna species considered except for albacore tuna.SST was the most important variable following year, month, and bottom depth for all the tuna species considered in this study.GAM analysis, modeling the CPUE of tuna species, revealed that catch year and catch month had biggest effect (Figure 4) on the fisheries success.
Although no scientific studies have been conducted on tuna habitat preferences in Tonga, areas with relatively high CPUE of tuna species could be attributed to oceanographic conditions associated the bathypelagic features within its EEZ.As previously stated, this central part of the EEZ partly envelops the Tonga Ridge, Tofua Arce Volcanic Front, the northern end of the Tonga Kermadec Arc, the westward region of the Lau Basin and the parallel north to south chains of volcanic seamounts along the Tonga Ridge.These oceanic features may influence environmental conditions, such as surface water temperature, nutrient, salt content, and mixed layer depth, which are habitats for pelagic species such as tuna.This may also be the reason for the persistent presence of the four tuna species throughout the year although the EEZ of Tonga is part of the WCPO and is generally characterized as oligotrophic [61], due to deep vertical mixing of water masses.Studies have demonstrated a strong correlation between tuna and shallow waters, particularly continental shelves and seamounts [62,63].These locations are widely recognized as prime habitats for large offshore fishes, primarily because of the substantial foraging advantages they offer [62] and possibly for reproductive and navigational benefits [64][65][66].
The temporal pattern in the standardized indices of abundance (Figure 7) and the partial effects plots suggested that during the time period of this study (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018), some of the stocks of the tuna population within the EEZ of Tonga experienced a gradual increase, whereas for other species it was mostly characterized by inter-annual variability.Some of the important regulatory changes and environmental dynamics within the time period of this study that could potentially be driving the trend in the CPUE of the four tuna populations within the EEZ include: (i) a fishing moratorium (2004-2011) on foreign fishing fleets [1], as part of the Tonga's long-term plan to expand its tuna fisheries, and (ii) the impact of El Niño Southern Oscillation (ENSO) events (such as those of 1997/98, 2003/04, 2009/10) known to affect tuna catchability through the spatial shift of tuna's preferred habitat away from normal fishing grounds in the Western Pacific [67,68].Strong seasonality was observed in the distribution and abundance of tuna populations.This could be attributed to seasonal variability, which is known to influence the migration patterns and habitats of tuna in the WCPO [5,13].Whilst Zainuddin et al., (2008) reported increased tuna abundance in cooler waters (April-September) in the Western North Pacific Ocean, Arrizabalaga et al. (2015) [20] reported higher abundance in warm waters due to tuna temperature preference for spawning.
Bathymetry also plays a role (Figure 5) in the central area of the EEZ of depths ranged between 2000 m and 3000 m (Figure S9), where the highest tuna capture occurred.As stated above, this is an area that partly envelops water of various depths due to its geographically topographic features.The eastern part of the EEZ covers the famous Tonga Trench, which runs south-north with depths of 5000 m or greater.
Our CPUE distributions (Figures S1-S4) show this area to have very low tuna encounters as compared to the central part.We attributed this habitat suitability distribution to two main reasons: (i) the productive waters in the deep areas are scattered over large vertical waters, whereas in smaller vertical waters, productive waters are compressed making them easy and quick to access for tuna, hence a higher tuna presence; (ii) since bathymetry is linked to both fleet operations and fish habitat; hence, it can act as a crude proxy for distance from ports.The distance from the main ports in Nuku'alofa (the capital of Tonga) to the high tuna occurrence areas (central areas of the EEZ) is in the range of 50 to 100 km, whereas the distance to the deepest areas is in the range of 150 to 400 km.Therefore, it could be that fishing vessels preferred the central area due to its shorter distance from ports, hence cutting back on operational costs.
The preferred ranges of the environmental factors, SST and Chl-a, are shown in Figure 4, but were found to be least important in modelling distribution, as measured by variable importance (values ≤ 0.1, Figure 5); much less important than spatial-temporal factors and bottom depth.Therefore, we tend to consider their influence over the tuna occurrence in the EEZ of Tonga to be limited (at least in the context of the present modelling framework).A study by Lan et al. (2017) [19] in the Tropical Pacific Ocean, suggests that areas with a higher SST and a Chl-a concentration of approximately 0.05-0.25 mgm −3 yield higher catch rates of yellowfin tuna.Setiawati et al. (2015) [10] indicated bigeye tuna habitat ranges for SST and Chl-a are 24.8 to 28.7 • C and 0.05 to 0.17 mgm −3 , respectively.
The results indicated by Zainuddin et al. (2006) that the distribution pattern of potential fishing zones for skipjack during the southeast monsoon were well-characterized by SST, ranging from 28.5 • C to 30.5 • C and with Chl-a ranging from 0.10 to 0.20 mgm −3 [55].In the Western Pacific Ocean, Zainuddin et al. (2008) report that the highest CPUE of albacore corresponded with areas of SST at 18.5-21.5• C and Chl-a of 0.2-0.4mgm −3 [69].In addition, studies have shown that areas with low catch encounters are found in the offshore waters due to the low Chl-a concentration present in these areas [69,70].Chl-a presence is an indicator of areas where small pelagic fishes concentrated for feeding [71].Hence, an adequate level of Chl-a concentration is required to enhance food availability and increase fish presence in the area.According to Atkinson et al. (2001) [72], Chl-a concentration of at least 0.2 mgm −3 is required (as compared to the Chl-a concentration in the EEZ of Tonga 0.03-0.1 mgm −3 ) to provide adequate levels of food to support a viable commerical fishery.Our results show that the highest Chl-a concentration in the archipelago waters of Tonga to be 0.1 mgm −3 , which is only half of that reported (0.2 mgm −3 ).Furthermore, there is no direct input of nutrient from land and coastal upwellings [61] to enhance primary productivity [73,74].
The spatiotemporal pattern in SST and Chl-a data (Figure S5 and Figure S6, respectively) show their monthly variability compared to maps of predicted suitable habitat for albacore (Figure 8), bigeye (Figure 9), skipjack (Figure 10), and for yellowfin (Figure 11).It was observed, as stated earlier, that the largest captures of the four tuna species are distributed in the central waters, 15 • S-24 • S, of the EEZ with temperatures above 22 • C and Chl-a concentration ≤ 0.001 mgm −3 .Our results show that albacore has the highest CPUE, followed by bigeye and yellowfin, and skipjack has the lowest as indicated by CPUE distributions (Figures S1-S4).The study by Arrizabalaga et al. (2015) [20] on global tuna catch found that albacore has the highest CPUE and is mostly caught between the latitudes of 20 • S and 40 • N, relatively high bigeye CPUE has been observed between 0 • and 40 • of each hemisphere, and yellowfin are found around the equator, while skipjack's highest CPUE are observed around 0 • -20 • of each hemisphere.The study by Hu et al. (2018) [75] found that bigeye tuna are caught more (more so than yellowfin, albacore, and skipjack) in equatorial waters farther from the coast and where the hypoxic layer is deeper.There was also indication of seasonal variation in their distributions, possibly corresponding to seasonal displacement of water temperature along low-high latitudes [76].Our results show they influence but do not seem to be the main determinants of the distribution and abundance of the tuna species in Tongan waters when compared, for example, to temporal measures like the variables month and year (Figure 5).

Conclusions
The study integrated catch data for four tuna species in Tonga's EEZ with satellitederived environmental data and employed distribution modeling and GAM analyses, revealing that high tuna catches align closely with predicted habitat maps from 2002 to 2018.The partial effect plots of our GAMs confirmed the non-linear relationships between the variables analyzed and tuna occurrence.Among the variables considered, SST exhibited the highest importance, followed by year, month, and bottom depth for all the tuna species in this study.The results of our study also suggest that the stocks of tuna populations within Tonga's EEZ showed a gradual increase for some species while others exhibited inter-annual variability.The key factors driving the CPUE trends for the four tuna populations could be attributed to a fishing moratorium on foreign fleets and the influence of El Niño-Southern Oscillation events altering the preferred habitat of tuna.Regarding spatial distribution, both the CPUE distribution and the predicted suitable tuna habitats in the study region indicate a preference for the tuna species to be concentrated in the central part of the EEZ, aligning with historically productive tuna fishing areas.The study revealed a seasonal peak in tuna occurrence between May and September, contrasting with lower occurrence rates from October to April, underscoring the consistent significance of spatiotemporal variables in shaping the habitat patterns of tuna and highlighting specific times and locations of abundance for these species.Overall, our model results predict the spatiotemporal pattern of the four tuna habitats in good agreement with our simple prediction maps and observation occurrence data.Future research needs to be conducted in the area in order to further examine the links and interactions between tuna habitats and other environmental and oceanographic factors not included in this study.

Diversity 2023 ,
15, x FOR PEER REVIEW 5 of 22 match the temporal scales of the environmental data.The distributions of the CPUEs of albacore, bigeye, skipjack, and yellowfin are shown in Figure 2, respectively.

Figure 2 .
Figure 2. Histograms of the distribution of CPUE (million tons per number of hooks) of albacore, bigeye, skipjack, and yellowfin for the duration of the study period of this study, 2002-2018, of the Tonga tuna fisheries within the EEZ, longitude 171.31 ° W-179.10 ° W.

Figure 2 .Figure 3 .
Figure 2. Histograms of the distribution of CPUE (million tons per number of hooks) of albacore, bigeye, skipjack, and yellowfin for the duration of the study period of this study, 2002-2018, of the Tonga tuna fisheries within the EEZ, longitude 171.31•W-179.10•W.

Figure 3 .
Figure 3. Histograms of the distribution of the environmental data, sea surface chlorophyll frontal gradients (chl_grad), depth, sea surface chlorophyll (log_chl), sea surface temperature and sea surface temperature frontal gradients (sst_grad) for the duration of the study period of this study, 2002-2018 of the NOAA repository [33], taken within the EEZ, longitude 171.31 • W-179.10 • W.

Figure 4 .
Figure 4.The partial effects for covariates chl_grad, depth, Chl-a, SST, sst_grad, month and year for albacore, bigeye, skipjack, and yellowfin, respectively, within the EEZ of Tonga, longitude 171.31 • W-179.10 • W. The partial effects of predictors for yellowfin tuna are presented in Figure 4.The partial effects of year indicates a gradual increase in CPUE from 2002 to 2018 with substantial inter- Diversity 2023, 15, x FOR PEER REVIEW 11 of 22

Figure 5 .
Figure 5. Index of relative importance (refer to Section 2.4) of predictors for the four tuna species in the EEZ of Tonga, longitude 171.31 ° W-179.10 ° W.Figure 5. Index of relative importance (refer to Section 2.4) of predictors for the four tuna species in the EEZ of Tonga, longitude 171.31•W-179.10•W.

Figure 5 .
Figure 5. Index of relative importance (refer to Section 2.4) of predictors for the four tuna species in the EEZ of Tonga, longitude 171.31 ° W-179.10 ° W.Figure 5. Index of relative importance (refer to Section 2.4) of predictors for the four tuna species in the EEZ of Tonga, longitude 171.31•W-179.10•W.

Figure 5 .
Figure 5. Index of relative importance (refer to Section 2.4) of predictors for the four tuna species in the EEZ of Tonga, longitude 171.31 ° W-179.10 ° W.

Figure 6 .Figure 6 .
Figure 6.Plots of posterior predictive distribution of the models for (A) albacore, (B) bigeye, (C) skipjack and (D) yellowfin in the EEZ of Tonga, longitude 171.31 ° W-179.10 ° W. Distribution of Figure 6.Plots of posterior predictive distribution of the models for (A) albacore, (B) bigeye, (C) skipjack and (D) yellowfin in the EEZ of Tonga, longitude 171.31 • W-179.10 • W. Distribution of the observed data is shown in black whereas the distributions from multiple draws of the model are shown in gray.Distribution of the observed data is shown in black whereas the distributions from multiple draws of the model are shown in gray.Posterior predictive check is way of checking if the model fits the data.The expectation is that if the model fits the data well, repeated generation of the data from fitted model should have the same distribution as the response variable.The presence of any systematic deviation from between the observed distribution and that posterior draws from the model indicates problem with the fitted model[48].In this case no issue was observed.

Figure 7 .
Figure 7. Time-series of standardized indices of CPUE for the four tuna species in the EEZ of Tonga, longitude 171.31 ° W-179.10 ° W.

Figure 7 .
Figure 7. Time-series of standardized indices of CPUE for the four tuna species in the EEZ of Tonga, longitude 171.31 • W-179.10 • W.

Table 2 .
Model summary for the best GAM model for the four species.The cumulative proportion of deviance explained (Cum.dev.explained) from the null model (a model with just an intercept term) to the full model is presented.