Skipjack tuna, Katsuwonus pelamis (Linnaeus, 1758), is a large pelagic species of significant economic importance in Indonesian fisheries, contributing substantially to national export revenues and serving as a major protein source for coastal communities (Firdaus 2019; Islamaji et al. 2020). Skipjack fisheries in Palabuhanratu Bay, West Java, Indonesia, are part of Fisheries Management Area (FMA) 573, which encompasses the southern waters of the Indian Ocean (Islamaji et al. 2020). The large pelagic fishery in the Palabuhanratu region remains highly promising for further development (Ihsan et al. 2017), positioning Palabuhanratu as a key base for small-scale fisheries that rely on tuna and skipjack resources. Notably, Palabuhanratu represents the largest bay-based fishing port complex in Java that directly faces the Indian Ocean, providing strategic access to productive oceanic fishing grounds and supporting intensive fishing activities in the region (Ihsan et al. 2017; Islamaji et al. 2020).

Skipjack tuna are commonly caught using various fishing gears, including payang (purse seine), troll lines, and longlines (Budiasih and Dewi 2015). Among these, troll lines are the most dominant and characteristic fishing gear employed by small-scale fishers in Palabuhanratu. This gear type is relatively simple and selective, resulting in minimal bycatch of non-target species (Chaliluddin et al. 2025). The success of trolling operations largely depends on the fishers’ ability to locate skipjack schools, whose distribution is dynamic in both spatial and temporal dimensions. In practice, fishers often rely on traditional knowledge and experiential “trial-and-error” methods, which can lead to low catch success and high fuel consumption (Nurdin et al. 2018).

Small-scale fisheries are particularly vulnerable to environmental variability due to their structural constraints, including limited capital investment, restricted access to navigational and forecasting technology, shorter operational range, and dependence on nearshore fishing grounds subject to highly dynamic coastal oceanographic conditions (Berkes et al. 2001; Islam and Chuenpagdee 2022). Unlike industrial fishing fleets that can relocate across broad oceanic areas in response to environmental shifts, small-scale troll-line fishers operating from ports such as Palabuhanratu are largely confined to coastal and neritic waters, where variations in oceanographic conditions can directly influence fish availability and catch success (Nurdin et al. 2018). In southern Java waters, environmental variability is associated with seasonal upwelling, riverine influences, and climate-driven processes linked to the Indian Ocean Dipole (IOD) and El Niño–Southern Oscillation (ENSO), all of which can affect marine productivity and fisheries performance (Susanto and Marra 2005; Iskandar 2014; Amalia et al. 2025). Consequently, improving the scientific understanding of how environmental variability shapes catch rates in these fisheries is not merely an academic exercise but represents a practical contribution to the adaptive capacity of small-scale fishing communities.

Skipjack tuna catch rates and availability are strongly influenced by oceanographic conditions. Factors such as sea surface temperature (SST) and chlorophyll-a (CHL) concentration are critical indicators of pelagic fish habitats, including skipjack (Sarianto 2018; Suhermat 2024). Other influential factors include salinity (Pratama et al. 2022) and sea level height (SSH) (Yati et al. 2024; Pratama et al. 2025a). Several studies have demonstrated that skipjack tend to associate with an optimal temperature range of 29.5–31.5°C and relatively high CHL concentrations due to enhanced food availability (Fauzan et al. 2018). Areas experiencing upwelling typically exhibit elevated nutrient concentrations that support primary productivity, thereby attracting fish aggregations (Susanto and Marra 2005).

However, marine ecosystems are influenced by interacting physical and biogeochemical processes that rarely operate in isolation. Consequently, analyses based on a limited number of environmental variables may overlook important ecological relationships associated with fish distribution and catch variability. Integrating multiple oceanographic indicators with fisheries data has proven useful for identifying potential fishing grounds and seasonal distribution patterns of pelagic species (Insanu et al. 2013; Pratama et al. 2022). For small-scale fishers in Palabuhanratu, the utilization of oceanographic information can potentially reduce search time, lower fuel costs, and improve catch efficiency (Nurani et al. 2021). Consequently, this approach not only enhances operational efficiency but also supports more informed and adaptive fisheries practices.

Although numerous studies have examined the influence of individual oceanographic parameters such as sea surface temperature and chlorophyll-a on skipjack tuna catches and catch per unit effort (CPUE), most previous research has relied on one or two variables, focused on industrial-scale fisheries, or covered relatively short temporal periods. In Palabuhanratu waters and surrounding areas along the southern coast of Java, existing studies have generally employed surface parameters such as sea surface temperature and chlorophyll-a (e.g., Nurdin et al. 2018; Pratama et al. 2025b; Suryapratama et al. 2025), while studies incorporating a broader range of variables remain limited to specific physical aspects such as bathymetry and salinity (Pratama et al. 2022). To date, research that simultaneously integrates physical and biogeochemical oceanographic parameters within a non-linear statistical framework remains limited. Consequently, the combined effects of physical–biogeochemical variability on skipjack tuna availability, particularly in small-scale troll line fisheries, are still not well understood.

This study addresses these gaps by integrating multi-year (2019–2023) physical and biogeochemical oceanographic data obtained from Copernicus Marine products with skipjack tuna CPUE derived from small-scale troll line fisheries operating from Palabuhanratu, Indonesia. Using a Generalized Additive Model (GAM), the study examines non-linear relationships between physical and biogeochemical oceanographic variables and skipjack tuna CPUE. The analysis focuses on identifying environmental conditions statistically associated with variations in catch rates, rather than defining skipjack tuna habitat suitability. Consequently, the results should be interpreted as statistical associations derived from fishery-dependent data that reflect fishery performance rather than causal ecological relationships.

Unlike most previous studies in the region, which have typically focused on one or two oceanographic variables or on industrial-scale fisheries, this study simultaneously evaluates the combined influence of multiple physical and biogeochemical indicators within a unified modelling framework using a multi-year dataset from a small-scale tropical fishery. By comparing the relative contribution of different environmental variable groups, the study provides new insights into how interacting oceanographic processes shape skipjack CPUE variability in tropical coastal waters. The findings contribute to the broader fisheries–oceanography literature and may support the development of more environmentally-informed fishing strategies and adaptive management approaches for small-scale fisheries facing increasing environmental variability.

Study area. This study was conducted to examine the relationship between physical and biogeochemical oceanographic variability and skipjack tuna catch per unit effort (CPUE) in the waters of Palabuhanratu, West Java, Indonesia. The analysis was conducted within the primary fishing grounds of the small-scale troll line fleet based at Palabuhanratu National Fishing Port (PPN Palabuhanratu), delineated as a fishing-ground polygon that served as the spatial framework for environmental data extraction and analysis (Fig. 1).

Fishing ground delineation and spatial data integration. Fishing location information was obtained from small-scale troll line fishery logbooks maintained at Palabuhanratu National Fishing Port (PPN Palabuhanratu), which recorded fishing positions in latitude and longitude coordinates. However, the availability of positional data was not consistent across all months during the study period. In addition, the coordinates recorded in the logbooks were only available at the degree and minute level without seconds, resulting in relatively limited spatial precision. Therefore, the fishing ground polygon was delineated using fishing locations documented during participatory field observations conducted in October–November 2022 and supported by fishermen’s GPS records as well as available positional information from the logbooks (Suppl. material 1: fig. S1). Although the field observations were conducted during a limited period, the resulting polygon was intended to represent the primary operational area of the small-scale troll line fleet based in Palabuhanratu rather than month-specific fishing locations. This assumption is supported by the available logbook coordinates from other years, which consistently fall within the same general area, and by the relatively localized nature of small-scale fisheries, whose operational range is commonly constrained by vessel size, fuel capacity, travel distance, and port-based fishing activities (Berkes et al. 2001; FAO 2015). Consequently, the polygon was used as a spatial framework for integrating fisheries and environmental datasets throughout the study period.

Subsequently, oceanographic parameters were extracted from all grid cells located within the fishing ground polygon. For each environmental variable, monthly values were calculated as the spatial average of all grid cells contained within the predefined fishing area. The resulting spatial aggregation generated a single representative monthly environmental value for each variable, which was subsequently paired with monthly CPUE data for GAM analysis. This approach enabled the consistent integration of fisheries and environmental datasets despite the incomplete availability of individual fishing-location records throughout the study period.

Data collection. Two primary datasets were used: fisheries production data and oceanographic data. Fisheries production data were obtained from official fishing logbooks maintained by Palabuhanratu National Fishing Port (PPN Palabuhanratu). These records constitute secondary fisheries data collected routinely by the port authority from small-scale troll-line vessels operating from Palabuhanratu. The logbooks contained information on monthly skipjack tuna landings [kg] and fishing effort, defined as the total number of fishing trips conducted each month. Monthly catch and effort data were subsequently aggregated for analysis. Catch per unit effort (CPUE) was calculated as the ratio of total monthly catch [kg] to the total number of fishing trips per month and was used as a relative index of catch rate.

Fishery-dependent CPUE may reflect not only fish abundance but also fishing behavior, spatial targeting strategies, vessel efficiency, and other operational factors (Maunder and Punt 2004). Accordingly, CPUE in this study is interpreted as an indicator of fishery performance and relative habitat association rather than an absolute measure of stock abundance. This distinction is particularly relevant in small-scale coastal fisheries where standardization of fishing effort is limited by data availability.

Oceanographic data were obtained from the Copernicus Marine Service (CMEMS), encompassing both physical and biogeochemical reanalysis products. Physical variables included sea surface temperature (SST), salinity (SAL), sea surface height (SSH), and current speed (CUR), while biogeochemical variables included chlorophyll-a concentration (CHL), pH, and silicate; nitrate, phosphate, and dissolved oxygen (DO) were additionally evaluated during the initial variable screening phase. All environmental variables used in this study were extracted from the surface layer of the respective Copernicus Marine products; therefore, no vertical integration or depth averaging was performed. Silicate was included not as a direct physiological driver of skipjack catch rates, but as a biogeochemical indicator of upwelling water mass intrusion. In coastal upwelling systems, dissolved silicate co-upwells with cold, nutrient-rich subsurface water and has been demonstrated as a reliable tracer of upwelling intensity (Li et al. 2011), while elevated silicate availability promotes diatom-dominated phytoplankton communities that underpin the short, energy-efficient food webs supporting skipjack prey fields (Goldenberg et al. 2024).

Physical variables were extracted from the Global Ocean Physics Reanalysis product GLORYS12V1 (https://doi.org/10.48670/moi-00021), which integrates satellite observations, in situ measurements, and numerical ocean circulation modelling through data assimilation to generate physically consistent estimates of ocean state variables (Lellouche et al. 2021). Biogeochemical variables were obtained from the Global Ocean Biogeochemistry Hindcast product (https://doi.org/10.48670/moi-00019), generated using the PISCES biogeochemical model within the NEMO modelling framework, forced by GLORYS12V1 ocean physics and ERA-Interim atmospheric reanalysis data. Unlike the physical reanalysis, the biogeochemical hindcast does not incorporate data assimilation, providing temporally consistent three-dimensional biogeochemical fields derived from model physics alone. These CMEMS products have been applied in regional oceanographic and fisheries studies in Indonesian and adjacent waters, including for skipjack tuna distribution modelling in the South Java Sea and Indian Ocean (Semedi et al. 2023; Pratama et al. 2025a), supporting their suitability for long-term analyses of environmental drivers. Detailed information on spatial resolution, temporal resolution, and product specifications is provided in Table 1.

Data preparation. Fisheries and oceanographic datasets were integrated into a single tabular dataset by aligning them to a common temporal resolution and spatial domain. Environmental variables from Copernicus Marine products with different native spatial resolutions, 0.083° (~9 km) for physical variables and 0.25° (~28 km) for biogeochemical variables, were processed independently. This approach does not constitute formal spatial resampling or reprojection between grids; instead, both datasets were independently aggregated to the same analytical unit, a single polygon-level monthly mean, regardless of their native resolution. It should be noted that the coarser 0.25° resolution of the biogeochemical product resulted in a substantially smaller number of grid cells contributing to the spatial averages within the fishing-ground polygon (12 grid cells) compared with the finer-resolution physical product (114 grid cells at 0.083° resolution), based on the centroid-in-polygon extraction procedure used throughout the analysis. Consequently, biogeochemical variables may be more sensitive to spatial averaging effects and may provide a less detailed representation of local environmental heterogeneity than the physical variables. This limitation is acknowledged as a source of uncertainty in the biogeochemical predictor estimates and should be considered when interpreting the results.

To improve data consistency, the analysis was restricted to troll line vessels within the 5–10 Gross Tonnage (GT) fleet segment operating from Palabuhanratu Fishing Port. This filtering approach was intended to reduce, although not completely eliminate, variability associated with differences in vessel size and fishing power among fleet segments. Nevertheless, some sources of uncertainty may remain because the selected fleet segment still encompasses vessels ranging from 5 to 10 GT, representing up to a twofold difference in vessel size. Vessels within this range can differ in operational characteristics, including engine capacity, crew size, fishing equipment, fishing strategy, and trip duration. These factors were not consistently available in the fisheries statistics and therefore could not be explicitly incorporated into the analysis.

Monthly environmental values were subsequently matched with monthly CPUE values. Data quality control included screening for extreme values and the removal of incomplete records. Prior to modeling, CPUE data were log-transformed [log(CPUE + 1)] to reduce right skewness and improve distributional properties for statistical analysis.

Before fitting the Generalized Additive Model (GAM), multicollinearity among the environmental predictor variables was carefully evaluated. Pairwise relationships were first examined using Pearson correlation analysis, with a threshold of |r| > 0.8 indicating strong collinearity (Sugiyono 2010). Variables exceeding this threshold were excluded based on ecological relevance and their contribution to model interpretability. Subsequently, the Variance Inflation Factor (VIF) was calculated for the remaining predictors to further assess collinearity structure and ensure model robustness, with a conservative threshold of VIF < 10 applied to minimize redundancy among explanatory variables and improve overall model stability and interpretability (Chen et al. 2025; Noor et al. 2025).

Generalized additive model (GAM) framework. The relationship between oceanographic variables and skipjack CPUE was analyzed using Generalized Additive Models (GAM). GAM was selected because of its ability to capture non-linear ecological responses commonly observed in fish–environment relationships (Zuur et al. 2009). Models were fitted using the mgcv package in R (Wood 2017), with Thin Plate Regression Splines (TPRS) applied as smoothing functions for all continuous predictors. GAM estimates smooth relationships between predictors and response variables through spline-based smoothing functions, allowing flexible representation of ecological responses without assuming a predefined functional form (Susilo and Wibawa 2017; Wood 2017). The Gaussian family was adopted because log-transformed CPUE values [log(CPUE + 1)] approximated a normal distribution. Model performance and variable contributions were evaluated using adjusted R², deviance explained, Akaike Information Criterion (AIC), estimated degrees of freedom (edf), and significance tests for smooth terms.

Model structure and selection. Three GAM schemes were developed to assess the relative contribution of different groups of oceanographic variables:

Scheme 1: Integrated physical and biogeochemical variables

Scheme 2: Physical variables only

Scheme 3: Biogeochemical variables only

This multi-scheme design allows comparison of the explanatory power of different oceanographic components in shaping skipjack CPUE variability (Table 2).

Model diagnostics and validation. Model assumptions were evaluated through residual diagnostics, including visual inspection of residual distributions, quantile–quantile (Q–Q) plots, and residual versus fitted value plots to assess normality and homoscedasticity (Lee 2025). Temporal autocorrelation in residuals was examined using autocorrelation function (ACF) plots.

Model performance was evaluated using the Akaike Information Criterion (AIC) and deviance explained (Wood 2017). AIC was used to compare model structures by balancing model fit and complexity, while deviance explained was used to quantify the proportion of variability in the response variable accounted for by the model (Susilo and Wibawa 2017). To assess model generalization, a k-fold cross-validation approach (k = 5) was applied, and prediction errors were compared across model schemes. This approach provides an evaluation of the model’s predictive performance on unseen data.

To further assess the contribution of individual predictors to overall model performance, a leave-one-variable-out (LOVO) sensitivity analysis was conducted for all environmental variables that remained statistically significant in the final GAM (salinity, SSH, chlorophyll-a, and pH). Separate models were fitted by sequentially removing one significant predictor at a time while retaining all remaining variables. Model performance was then compared with the full model using AIC, deviance explained, and adjusted R². Changes in these metrics relative to the full model were used to evaluate the relative contribution of each predictor and to assess whether model performance was disproportionately influenced by any single environmental variable.

By combining these complementary metrics, the most suitable model was identified as the one that best balances goodness-of-fit, explanatory power, and predictive performance in explaining the relationship between oceanographic variables and skipjack CPUE.

Visualization and interpretation. Partial response curves were generated to visualize the non-linear effects of each oceanographic variable on skipjack CPUE. The x-axis represents the range of predictor variables, while the y-axis indicates their partial effects on CPUE. Environmental ranges associated with higher CPUE were identified based on peak response values, which represent conditions associated with higher skipjack CPUE.

Methodological limitations. This study has several limitations that should be considered when interpreting the results. First, environmental variables were obtained through spatial aggregation by averaging all grid cells located within the fishing ground polygon for each month of observation. This approach was adopted to match the temporal resolution of the environmental data with the available CPUE data and to represent the average environmental conditions across the fleet’s primary fishing area. However, such spatial aggregation may reduce the influence of fine-scale oceanographic features, including thermal fronts, eddies, upwelling cells, and other habitat structures that may affect skipjack tuna distribution at finer spatial scales. Consequently, the environmental variables used in the GAM represent average environmental conditions across the main fishing grounds rather than the conditions directly experienced by fish at individual fishing locations.

Second, the use of monthly aggregated data may not fully capture short-term environmental variability (daily to weekly) that could be relevant to fishing operations. Temporal aggregation may therefore obscure transient oceanographic events that potentially influence fish availability and catch rates.

Furthermore, the CPUE data used in this study were derived from fishery-dependent sources and may therefore reflect not only environmental influences but also fishing behavior and operational strategies. Although the analysis was restricted to troll-line vessels within the 5–10 GT fleet segment to reduce heterogeneity among vessels, variation in engine capacity, crew size, gear configuration, fishing strategy, and trip duration could not be explicitly standardized because such information was not consistently available. Therefore, the results should be interpreted as broad-scale environmental associations related to fishery performance rather than as direct indicators of habitat suitability or stock abundance.

Temporal variability of skipjack tuna CPUE. Descriptive analysis of catch per unit effort (CPUE) and number of trips (effort) data for skipjack tuna, Katsuwonus pelamis, during the 2019–2023 period revealed substantial temporal variability (Fig. 2). Monthly CPUE values ranged from approximately 67 kg per trip in July 2019 to a peak of 778 kg per trip in July 2022. Higher CPUE values were generally observed during 2022–2023 compared with the earlier years of the study period, although considerable month-to-month fluctuations occurred throughout the time series. These results indicate that skipjack CPUE varied markedly across both months and years during the study period.

The temporal pattern of CPUE showed recurring fluctuations throughout the study period, with relatively higher catch rates generally occurring during the mid-year months (June–September). Particularly high CPUE values were recorded during July and August 2022, which represented the highest catches observed during the study period. In contrast, several months during 2019 were characterized by comparatively lower CPUE values. Overall, the observed patterns suggest the presence of both seasonal and interannual variability in skipjack catch rates in Palabuhanratu waters.

Distribution of oceanographic parameters. The distribution of physical oceanographic parameters, including sea surface temperature (SST), salinity, sea surface height (SSH), and current speed, as well as biogeochemical parameters, including chlorophyll-a concentration, pH, silicate, nitrate, phosphate, and dissolved oxygen (DO), is presented using boxplots (Fig. 3).

Descriptive analysis of oceanographic parameters revealed substantial variation among variables. SST ranged from 25.49 to 29.94°C, with a median of approximately 28.7°C, indicating warm tropical water conditions. Salinity was relatively stable, ranging from 33.10 to 34.47‰¹ with a median of 33.7‰. SSH varied between 0.28 and 0.67 m (median around 0.55 m), while surface current speeds ranged from 0.10 to 0.42 m · s^–1 with a median of 0.23 m · s^–1, suggesting predominantly calm to moderate currents.

CHL concentration was generally low, with a median of 0.25 mg · m^–3, although several high values (>0.9 mg · m^–3) appeared as outliers. Similarly, silicate ranged from 1.15 to 3.9 mmol · m^–3 with a median of 2.6 mmol · m^–3. Other nutrient parameters, such as nitrate and phosphate, exhibited very low concentrations (median close to 0), yet occasional high outliers indicated episodic nutrient input events.

DO ranged from 199 to 212 mmol · m^–3 with a median of 201 mmol · m^–3, while pH remained relatively stable within a narrow range of 7.99–8.05 (median = 8.03). In contrast to the relatively limited variability observed for DO and pH, several nutrient-related variables exhibited broader distributions. The presence of extreme values in CHL, nitrate, and phosphate suggests periods of enhanced biological productivity and variability in nutrient availability, which may reflect seasonal changes in oceanographic conditions within the study area.

Multicollinearity test. The results of the multicollinearity analysis (Fig. 4) indicate that several oceanographic variables exhibited very high correlations (r > 0.8), potentially leading to multicollinearity in the model. Based on the heatmap, strong correlations were observed between variable pairs such as CHL–Nitrate (r = 0.91), CHL–Phosphate (r = 0.95), Nitrate–Phosphate (r = 0.96), DO–CHL (r = 0.91), DO–Phosphate (r = 0.94), and DO–Nitrate (r = 0.85). These findings indicate that CHL, Nitrate, Phosphate, and DO belong to a highly collinear group and therefore must be excluded from the model to avoid redundancy and improve model robustness.

To avoid bias and maintain ecological interpretability, variable selection was based on both statistical screening and ecological relevance. Among the highly collinear biogeochemical variables identified in the correlation analysis, CHL was retained as the representative predictor to minimize redundancy in the final model. In addition to exhibiting greater temporal variability than the nutrient variables, CHL directly reflects phytoplankton biomass and marine productivity, which are commonly associated with pelagic fish distribution and catch variability (Usman et al. 2017), and could be independently evaluated against MODIS OceanColor observations. Consequently, nitrate, phosphate, and dissolved oxygen (DO) were excluded from subsequent modelling. Therefore, the final model focused on variables that were both statistically defensible and ecologically meaningful while minimizing redundancy among predictors and improving model interpretability.

Multicollinearity diagnostics showed that all predictor variables had acceptable VIF values (<10) (Table 3). Moderate collinearity observed among SST, CHL, and pH reflects natural coupling between physical and biogeochemical processes in coastal environments. All variables were therefore retained based on their ecological relevance and acceptable statistical properties.

GAM model results. The Generalized Additive Model (GAM) constructed after eliminating nitrate, phosphate, and DO, produced a simpler structure with seven oceanographic variables: salinity, SST, SSH, current speed, CHL, pH, and silicate. The removal of the three highly collinear variables improved model stability and estimation reliability. The final model explained 70.7% of deviance, with an adjusted R² of 0.623, indicating strong model performance (Table 4).

Based on statistical significance, salinity emerged as the most dominant factor, exhibiting a strong non-linear relationship (edf = 1.584; F = 19.617; P < 0.05). In addition, CHL (edf = 1.000; F = 8.217; P < 0.01), pH (edf = 1.000; F = 5.419; P < 0.05), and SSH (edf = 3.200; F = 3.324; P < 0.05) also showed significant effects on CPUE. Meanwhile, SST showed a marginally significant effect (edf = 2.692; F = 2.433; P = 0.0679). Current speed (P = 0.109) and silicate (P = 0.215) did not show significant effects in the model.

GAM model evaluation. Model performance was evaluated using complementary approaches to assess both model fit and generalization ability. The Akaike Information Criterion (AIC) was used to compare model structures, while deviance explained was used to quantify the proportion of variability captured by the model (Wood 2017). Lower AIC values indicate a better balance between model fit and complexity.

A comparison of AIC values among the three GAM schemes (Table 5) shows that the model integrating both physical and biogeochemical variables produced the lowest AIC (743.26), outperforming the physical-only (762.64) and biogeochemical-only (778.73) models. This indicates that combining multiple environmental components provides a more comprehensive explanation of skipjack CPUE variability.

In addition to AIC-based evaluation, model generalization was assessed using k-fold cross-validation (k = 5). The results showed variable prediction errors across folds, with root mean square error (RMSE) values ranging from 0.356 to 0.749 and a mean RMSE of 0.485 ± 0.137. The corresponding mean absolute error (MAE) was 0.369 ± 0.101, indicating moderate prediction error on the log-transformed CPUE scale. This variability likely reflects the heterogeneous and dynamic nature of fisheries data, while the overall prediction error indicates that the model maintains moderate predictive performance without strong evidence of overfitting.

A leave-one-variable-out (LOVO) sensitivity analysis further supported the robustness of the final model (Suppl. material 1: table SS4). Excluding salinity produced the largest deterioration in model performance, increasing AIC by 24.01 units and reducing deviance explained from 70.70% to 54.63%. Excluding SSH, chlorophyll-a, and pH also resulted in higher AIC values and lower explanatory power relative to the full model, although the magnitude of change was substantially smaller (ΔAIC = 4.61–5.35). These results indicate that salinity contributed the greatest share of explanatory information, while all significant predictors provided complementary information that improved model performance. Collectively, the sensitivity analysis suggests that the final GAM was not driven by a single predictor variable but rather by the combined influence of multiple environmental factors.

Residual diagnostics further confirmed that model assumptions were reasonably satisfied. Visual inspection of residual distributions and Q–Q plots indicated no substantial deviation from normality, while residual versus fitted plots did not reveal clear heteroscedastic patterns (Lee 2025). Autocorrelation analysis using ACF plots showed no strong temporal dependence in residuals, suggesting that the independence assumption was adequately met.

Detailed results of model validation and diagnostic analyses, including k-fold cross-validation metrics, leave-one-variable-out sensitivity analysis, Q–Q plots, and autocorrelation function (ACF) plots, are provided in the Suppl. material 1: figs S2–S4.

Response curve. The GAM results indicate that skipjack tuna CPUE exhibited different responses to individual oceanographic parameters, suggesting that CPUE variability is associated with a combination of multiple environmental factors. Based on model evaluation results (Table 5), the model integrating both physical and biogeochemical oceanographic parameters (Scheme 1) produced the lowest AIC value and was therefore selected as the final model. Consequently, the response curves (Fig. 5) represent the partial effects of each environmental variable derived from the final GAM after accounting for the influence of the other predictor variables included in the model.

Salinity exhibited a negative relationship with CPUE, with CPUE tending to decrease as salinity increased. SST also showed a generally negative relationship with CPUE. SSH displayed a negative trend, although uncertainty increased toward the extremes of the observed data range. Current speed exhibited a non-linear response, with relatively higher CPUE values associated with low to moderate current velocities and lower CPUE values at higher current speeds. Among the biogeochemical variables, CHL showed a positive relationship with CPUE. In contrast, pH exhibited a negative relationship with CPUE. Silicate also displayed a generally negative trend in relation to CPUE. Overall, the response curves suggest that skipjack tuna CPUE is associated with a combination of physical and biogeochemical oceanographic parameters, each exhibiting distinct response patterns.

Skipjack tuna, Katsuwonus pelamis, CPUE in Palabuhanratu waters during 2019–2023 exhibited distinct temporal fluctuations, with higher catch rates recorded during 2022–2023 relative to preceding years and a peak of 778 kg per trip in mid-2022 (Fig. 2). Concurrent temporal variability in several oceanographic variables (Fig. 3; Suppl. material 1: table SS1) suggests that shifts in environmental conditions contributed to the observed patterns in catch rates.

The study period encompassed contrasting phases of the Indian Ocean Dipole (IOD). In 2023, a positive IOD phase developed from August to December; positive IOD conditions are generally associated with enhanced upwelling intensity in the eastern tropical Indian Ocean, elevating nutrient availability, reducing SST, and increasing CHL concentrations (Fadholi 2013; Iskandar 2014). In contrast, 2022 was dominated by a negative IOD phase, which is typically associated with higher SST, weaker upwelling, and reduced nutrient supply to surface waters (Martono 2017; Nurani et al. 2021). Notably, relatively high CPUE values were observed during both IOD phases, indicating that catch variability at Palabuhanratu cannot be attributed to IOD forcing alone. Instead, CPUE variability is likely shaped by the interaction of multiple oceanographic processes at regional and local scales, including possible modulation by ENSO (Amalia et al. 2025). It should be noted that the IOD index was not included as a predictor in the GAM; its discussion here provides regional climatic context only.

The GAM results (Table 4) confirmed that skipjack CPUE was associated with a combination of physical and biogeochemical variables. The integrated model incorporating both variable groups outperformed models based on either group alone, producing the lowest AIC value (Table 5), and demonstrated the capacity of the GAM approach to capture non-linear relationships characteristic of ecological processes in dynamic marine systems. A leave-one-variable-out sensitivity analysis (Suppl. material 1: table SS4) further confirmed that no single predictor dominated model performance: excluding salinity produced the largest deterioration, increasing AIC by 24.01 units and reducing deviance explained from 70.7% to 54.6%, while excluding SSH, CHL, and pH also degraded model fit relative to the full model. Together, these results confirm that both physical and biogeochemical variables contributed unique and complementary explanatory information to the integrated GAM.

Salinity emerged as the most influential predictor (P < 0.001), with CPUE tending to decrease as salinity increased across the observed range (Fig. 5). This pattern is consistent with salinity’s role as an indicator of water-mass origin and mixing dynamics. Lower surface salinity in Palabuhanratu waters may reflect terrestrial freshwater inputs from river systems such as the Cimandiri River, or upwelling-driven advection of coastal water masses, both of which are associated with enhanced nutrient availability and productivity. This association between lower salinity and higher catch rates is consistent with previous studies identifying salinity as an important environmental correlate of pelagic fish occurrence in regional waters (Laevastu and Hayes 1981; Baskoro et al. 2011; Ma’mun et al. 2019; Pratama et al. 2022).

Sea surface temperature (SST) exhibited a negative association with CPUE, with higher catch rates tending to occur under cooler surface conditions (Fig. 5), though this relationship was only marginally significant (P = 0.068). In the southern Java upwelling system, cooler SST is commonly associated with the advection of nutrient-rich subsurface water and the development of frontal structures that enhance prey aggregation (Nurdin et al. 2018; Prasetya et al. 2022). Rather than representing a direct physiological constraint, SST in this context likely acts as an integrative indicator of broader oceanographic productivity processes that structure foraging conditions for skipjack tuna.

Lower SSH values were associated with higher CPUE, consistent with the role of negative SSH anomalies as indicators of upwelling, convergence zones, and mesoscale features such as fronts and eddies that concentrate prey organisms and attract pelagic predators (Semedi et al. 2023; Pratama et al. 2025a). Although current speed was not statistically significant, the response curve suggests that moderate current speeds may support favorable foraging conditions by facilitating prey aggregation without imposing excessive energetic demands on skipjack tuna (Jalil 2013).

From a biogeochemical perspective, CHL showed a positive association with CPUE, with higher catch rates at elevated concentrations (Fig. 5). CHL is a widely used proxy for phytoplankton biomass and primary productivity, and its positive relationship with pelagic predator occurrence is well-documented in tropical and subtropical marine systems (Parsons et al. 1984; Itoh et al. 2011; Pratama et al. 2025a). In the neritic and coastal waters of the fishing-ground polygon, CHL variability likely integrates signals from both upwelling-driven nutrient enrichment and riverine inputs from river systems surrounding Palabuhanratu Bay, particularly the Cimandiri River. The high coefficient of variation for CHL (CV = 110.48%; Suppl. material 1: table SS1) confirms that this variable captured substantial and ecologically meaningful temporal variability throughout the study period, making it a particularly informative predictor of foraging conditions for skipjack tuna.

pH was statistically significant in the model; however, its very narrow observed range (7.99–8.05; CV = 0.17%) substantially limits confidence in the ecological interpretation of its response. The significant association likely reflects underlying biogeochemical co-processes rather than a direct physiological effect on skipjack tuna. In productive neritic environments, slight variations in pH are associated with photosynthesis–respiration dynamics and freshwater–seawater mixing (Heuer and Grosell 2014), both of which intensify in coastal systems influenced by terrestrial inputs, reinforcing pH’s interpretation as a broad-scale proxy for biogeochemical activity rather than a primary ecological driver. Given the coarser spatial resolution (0.25°) and absence of data assimilation in the biogeochemical hindcast product, its ecological interpretation warrants corresponding caution.

Silicate did not exhibit a statistically significant relationship with CPUE (P = 0.215, edf = 1.000). Silicate was included as a biogeochemical indicator of upwelling water-mass intrusion, given that it co-upwells with cold, nutrient-rich subsurface water and serves as a recognised tracer of upwelling intensity in coastal systems (Li et al. 2011). In upwelling environments, elevated silicate supports diatom-dominated phytoplankton communities that underpin energy-efficient food webs for small pelagic prey of skipjack tuna (Turner et al. 1998; Goldenberg et al. 2024). The absence of a significant response likely reflects that CHL already captures the integrated productivity signal of upwelling in this system, rendering the independent contribution of silicate statistically indistinguishable at monthly resolution. This finding highlights the value of evaluating biogeochemical variables beyond CHL in future studies, particularly with higher temporal resolution or direct prey biomass measurements.

The contrasting levels of temporal variability among predictors provide important context for interpreting the model results. While pH (CV = 0.17%) and salinity (CV = 1.02%) exhibited limited relative variability, other variables showed substantially greater temporal variation: CHL (CV = 110.48%), current speed (CV = 34.73%), silicate (CV = 21.11%), and SSH (CV = 19.96%) (Suppl. material 1: table SS1). These differences indicate that environmentally meaningful variability was present throughout the study period despite the narrow ranges observed for some parameters and reinforce the need to interpret low-CV variables such as pH with greater caution. The finding that the integrated model outperformed single-domain alternatives (Table 5) is consistent with this picture, as it reflects the complementary environmental signals contributed by variables with contrasting variability structures.

Oceanographic variables were extracted and averaged within a fishing-ground polygon constructed from available fishing-location records (Suppl. material 1: fig. S1), as individual trip-level coordinates were not consistently available across the full study period. Although this approach enabled the systematic integration of fisheries and gridded oceanographic datasets, spatial averaging is known to smooth environmental gradients and reduce representation of localized features such as thermal fronts and mesoscale eddies (Wiens 1989; Fotheringham and Wong 1991). Accordingly, the relationships identified here should be interpreted as associations between CPUE and average environmental conditions within the main fishing grounds, rather than fine-scale habitat relationships at individual fishing locations.

The limitations of global reanalysis and hindcast products also warrant consideration. Copernicus Marine products are designed for regional to global applications and may not fully resolve environmental variability in complex coastal upwelling environments, particularly for temperature and salinity (de Souza et al. 2020; Castillo-Trujillo et al. 2023), with local-scale variability further dampened by data assimilation and modelling procedures. Nevertheless, several lines of evidence support the suitability of these products for the present analysis. They have been widely applied in fisheries oceanography research in Indonesian waters, including GAM-based studies of skipjack tuna in southern Java (Semedi et al. 2023), and the GLORYS12 reanalysis underpinning several Copernicus physical products was specifically designed for physically consistent regional ocean applications (Lellouche et al. 2021). Independent validation against RAMA buoy in-situ observations in Indonesian waters demonstrated strong agreement for physical variables, with correlation coefficients of 0.988, 0.960, and 0.934 for SST, SSS, and SSH, respectively (Nurzeha et al. 2025). Additional cross-comparison against MODIS OceanColor observations during 2019–2023, conducted in the present study, confirmed strong temporal agreement for SST (Pearson r = 0.938, Spearman ρ = 0.913; P < 0.001) and CHL (Pearson r = 0.886, Spearman ρ = 0.725; P < 0.001), with RMSE values of 0.897°C and 1.261 mg · m^–3, respectively (Suppl. material 1: table SS2). Although absolute differences were evident, particularly for CHL, the consistently high Pearson and Spearman correlations indicate that Copernicus products reliably reproduced the timing, direction, and relative magnitude of monthly environmental fluctuations throughout the study period. This distinction is important because the objective of the present analysis was not to estimate absolute environmental conditions at individual fishing locations, but rather to evaluate statistical associations between temporal environmental variability and CPUE at the monthly scale. Consequently, the ability of the Copernicus products to capture coherent temporal patterns is of greater relevance to the GAM framework than exact agreement in absolute values. The larger RMSE observed for CHL likely reflects differences in retrieval methodology, spatial resolution, and data-processing approaches between satellite-derived and model-based products and therefore supports interpreting CHL as a broad-scale indicator of ecosystem productivity rather than a precise point measurement. For pH and silicate, which were not independently cross-validated and are derived from a coarser-resolution (0.25°) biogeochemical hindcast without data assimilation, a greater degree of uncertainty should be assumed, consistent with the cautious interpretation adopted throughout this study.

Model evaluation results further support the robustness of the GAM framework. The integrated model explained approximately 70%–71% of the deviance, achieved an adjusted R² exceeding 0.60, and produced the lowest AIC among the candidate models. k-fold cross-validation (k = 5) yielded moderate and reasonably consistent prediction errors across folds (mean RMSE = 0.485 ± 0.137; mean MAE = 0.369 ± 0.101), reflecting the inherent heterogeneity of fishery-dependent CPUE data (Maunder and Punt 2004). The consistency of prediction performance across folds, combined with high deviance explained, suggests that the model captured recurrent environmental associations rather than overfitting to random fluctuations in the dataset.

It is important to emphasize that the relationships identified in this study represent statistical associations rather than evidence of direct ecological causation. CPUE as a fishery-dependent index may be influenced by fishing behavior, spatial targeting, and operational factors beyond environmental conditions (Maunder and Punt 2004). Monthly-averaged oceanographic data may likewise fail to capture short-term environmental variability relevant to individual fishing operations. The identified environmental ranges should therefore be understood as conditions statistically associated with relatively higher catch rates under the analytical framework adopted, not as definitive characterizations of skipjack tuna habitat preferences.

In summary, skipjack CPUE variability in Palabuhanratu waters during 2019–2023 was associated with the combined influence of multiple physical and biogeochemical processes, with salinity, SSH, CHL, and pH emerging as the most informative predictors in the integrated GAM. The superior performance of the integrated model relative to single-domain alternatives underscores the importance of simultaneously considering physical and biogeochemical environmental dimensions in fisheries–oceanography research. These findings reflect the complexity of tropical coastal marine ecosystems and suggest that future studies would benefit from higher-resolution environmental observations, spatially explicit fishery data, and fishery-independent abundance indices capable of resolving ecological processes at finer spatial and temporal scales.

This study demonstrates that variations in skipjack tuna CPUE in Palabuhanratu waters are associated with a combination of physical and biogeochemical oceanographic conditions, as captured by a Generalized Additive Model (GAM). The best-performing model explained 70.7% of the CPUE deviance, indicating good performance in representing catch-rate variability. Higher CPUE values were generally associated with lower sea surface temperature, lower salinity, and higher chlorophyll-a concentrations within the observed range of environmental conditions. In addition, pH was identified as a significant biogeochemical predictor, although its ecological interpretation warrants caution due to its very narrow observed range and the greater uncertainty associated with the biogeochemical hindcast product. Furthermore, model selection results showed that the model integrating both physical and biogeochemical oceanographic variables performed better than models based on either variable group alone. These findings provide evidence that catch-rate variability in small-scale skipjack fisheries in tropical waters is not associated with a single environmental factor but is instead influenced by the interaction of multiple oceanographic processes operating simultaneously. Consequently, this study contributes to strengthening our understanding of the links between marine environmental variability and the dynamics of tropical pelagic fisheries.

Nevertheless, it should be emphasized that the results represent statistical associations derived from fishery-dependent data and should not be interpreted as direct indicators of habitat suitability or as evidence of causal ecological relationships. Future studies integrating higher spatial and temporal resolution datasets, together with independent fishery observations, are needed to improve understanding of habitat dynamics and the ecological mechanisms underlying variability in skipjack tuna catch rates.