CPUE standardization for southern bluefin tuna (Thunnus maccoyii) in the Korean tuna longline fishery, accounting for spatiotemporal variation in targeting through data exploration and clustering

Data

Set by set catch and effort data were compiled by the Korean National Institute of Fisheries Science (NIFS). Data were selected with the criterion that when a vessel reported the capture of at least one SBT in a month, all effort for the vessel-month was included.

The fields reported in the operational set by set data were vessel identifier (call sign), fishing location to 1° cell of latitude and longitude, date, effort (number of hooks and floats), and catch in numbers of southern bluefin tuna (SBT), bigeye (BET), yellowfin (YFT), albacore (ALB), skipjack (SKJ), swordfish (SWO), black marlin (BLM), blue marlin (BUM), striped marlin (MLS), sailfish (SFA), sharks (SHA), and other species (OTH).

Data used in this study were from 1996 to 2018. Data prior to 1996 were not available due to insufficient information for CPUE standardization. Dates were converted to months and quarters. Since Korean longliners set at night or at dawn, moon phase was used to calculate the relative lunar illumination for each date, using the R package lunar (Lazaridis, 2014). Spatial positions were classified into 5° cells, and CCSBT statistical areas (CCSBT, 2019b). The numbers of hooks between floats (HBF) were calculated by dividing hooks by floats and rounding to the nearest whole number.

For CPUE standardization, data were cleaned by removing sets in which there were fewer than 1,000 hooks or more than 5,000 hooks. Korean tuna longline vessels fishing for SBT in the CCSBT convention area have individual annual quotas. The fishing season is from April to March of the following year. They have mainly operated in two locations to the south of 35°S, either between 10°W and 50°E (within CCSBT statistical area 9) or between 90°E and 120°E (within CCSBT statistical area 8) (Fig. 2). Effort has focused on the western area (statistical area 9) from March to September/October and shifted to the eastern area (statistical area 8) from July/August until December (Fig. 3). For that reason, we defined two separate core SBT fishing areas: within statistical area 9 from March to October and statistical area 8 from July to December.

Data exploration

Data were plotted to explore the spatial and temporal distributions of effort, and patterns in operational characteristics such as the hooks per set and HBF. Operational characteristics were compared with catch rates to identify possible gear-based criteria for targeting. We examined patterns through time and among major species in the nominal catch rates by year-quarter and area and compared them with patterns in the proportions of sets with no catch of each species. We plotted maps of the species composition through time, to identify possible spatial and temporal variation in fishing behaviour or population composition.

To further explore changes in the fishery and identify periods of change, we plotted the participation of vessels in the fleet, sorted first by the start date and then by the end date of participation in the fishery. To explore changes in effort distribution and concentration through time, we plotted the numbers of 5° × 5° and 1° × 1° cells fished and the average number of operations per fished cell for each year and area.

Target change

Target change can be a significant problem for CPUE standardization since it can bias CPUE trends. Analyses were carried out using two alternative approaches to differentiate targeting practices.

GLM analyses

SBT CPUE was standardized using the set by set data and generalized linear models (GLMs) in Microsoft R Open 3.3.2 (R Core Team, 2016), and the methods generally followed the approaches used by Hoyle & Okamoto (2011) and Hoyle et al. (2015). Analyses were conducted separately for each of the two core areas, and for each of the two target change methods.

Data were prepared by selecting data for vessels that had made at least 100 sets, for years in which there had been at least 100 sets, and for 5° cells in which there had been at least 200 sets. Categories with too few sets provide estimates with high uncertainty and low reliability, so this approach removes a few cells, vessels, and time periods that lack much usable information.

The CPUE standardization was carried out using generalized linear models (GLMs) with both a lognormal constant model and a delta lognormal approach. The lognormal constant GLM was used to summarize the effects of covariates on the index (via the package influ, Bentley et al., 2011) across the whole dataset, but was not used for inferences about the abundance trend, since this approach has been superseded by methods that model zeroes more directly (Maunder & Punt, 2004). The preferred abundance indices were obtained using the delta lognormal approach.

Covariates for all models were specified as $\mathrm{year}+\mathrm{vessid}+\mathrm{latlong}+\lambda \left(\mathrm{hooks}\right)+g\left(\mathrm{month}\right)+h\left(\mathrm{moon}\right)$. The functions λ, g and h were cubic splines with 5, 3, and 4 degrees of freedom, respectively, with sufficient flexibility to explain variability while avoiding overfitting. Higher order terms are inadvisable for conventional polynomials but perform relatively well with regression splines (Harrell, 2001). The number of hooks (hooks) was included in the model to allow for possible hook saturation and potential targeting changes associated with hooks per set. The variable moon was the lunar illumination on the date of the set, which was included to find out whether SBT catch rate is related to moon phase. The variables year, vessid, and latlong (5° latitude-longitude cell) were fitted as categorical variables. For the clustering-based approach, the cluster was also included as a categorical variable.

The models did not include HBF directly. The data selection method had already addressed HBF by including only a narrow range of HBF values in the range 9–12. The cluster analysis method addresses targeting independently of HBF, and in any case, less than 1% of sets included HBF outside the 9–12 range.

The following lognormal model was used. ${\displaystyle ln\left(\mathrm{CPUE}+k\right)\sim \mathrm{covariates}+ϵ}$

The units of the input CPUE is catch-in-number of SBT per 100 hooks, and the constant k, added to allow for modelling sets with zero catches of the species of interest, is 10% of the mean CPUE for all sets (Campbell, 2004).

The delta lognormal approach (Lo, Jacobson & Squire, 1992; Maunder & Punt, 2004) used a binomial distribution for the probability w of catch rate being zero and a probability distribution f (y), where y was log (catch/hooks per set), for non-zero (positive) catch rates.

${\displaystyle Pr\left(Y=y\right)=\left\{\begin{array}{cc}w\hfill & y=0,\hfill \\ \left(1-w\right)f\left(y\right)\hfill & \mathrm{otherwise}\hfill \end{array}\right.}$ ${\displaystyle g\left(w\right)=\left(\mathrm{CPUE}=0\right)\sim \mathrm{covariates}+ϵ}$ ${\displaystyle f\left(y\right)=\mathrm{CPUE}\sim \mathrm{covariates}+ϵ}$

where g is the logistic function.

Data in the models were ‘area-weighted’, with the statistical weights of the sets adjusted so that the total weight per year in each 5° cell would sum to 1. This method was based on the approach identified using simulation by Punsly (1987) and Campbell (2004), that for set j in area i and year t with hooks h_ijt, the weighting function that gave the least average bias was: $w_{ijt}=\frac{\log \left(h_{ijt}+1\right)}{{\sum }_{j=1}^{n}log\left(h_{ijt}+1\right)}$. Given the relatively low variation in number of hooks between sets in a stratum, we simplified this to $w_{ijt}=\frac{h_{ijt}}{{\sum }_{j=1}^n{h}_{ijt}}$.

The models had no interactions between year effects and other covariates, so relative annual expected responses for the lognormal constant index and the lognormal component of the delta lognormal index were invariant with different values of other covariates. We generated an index for each of these models by predicting the response for each year with covariate values held constant.

For the delta model, however, the annual trend is affected by the values chosen for the other covariates, which are held constant when predicting annual catch rates. Choosing covariate values that give a higher rate of nonzero catch will reduce the variability among years in the delta index, and hence in combined index. To avoid subjectivity, the constant used to predict from the delta regression was adjusted so that the mean of the annual proportions of positive catches was the same in the predictions as in the observed data.

The combined index estimated for each year was the product of the year effects for the two model components, $\left(1-w\right)\times E\left(y|y\ne 0\right)$. This index was normalized to average 1, so the final index represents relative catch rate.

Model fits were examined by plotting the residual densities and using Q–Q plots.

Data exploration

Almost all Korean tuna longline vessels fishing for SBT used between 9 and 12 HBF (Fig. 4), with the majority of HBF outside this range coming from north of 35°S, outside the main SBT targeting area. The number of hooks per set has been relatively consistent since 1996, averaging a little over 3,000.

Mean nominal catch rates in statistical areas 8 and 9 were higher for SBT than for other species until the mid-2000s (Fig. S1). After this time in the areas 8 and 9 the SBT catch rates decreased and other species, particularly albacore tuna (ALB), increased. However, in the most recent years the SBT catch rates were again higher than other species. Similarly, the proportion of sets reported with zero SBT catches was low through most of the time series in the areas 8 and 9, but increased from 2004 to 2010 in area 9 and in some years during the early 2010s in area 8 (Fig. S2).

Statistical areas 2 and 14 in the Indian Ocean are at temperate latitudes between 20°S and 35°S. Highest catch rates were for YFT and more recently ALB in area 14, and BET and ALB in area 2 (Fig. S1). Since the mid-2000s ALB catch rates have increased markedly and particularly in the area 2, suggesting a trend towards targeting this species. Catch rates of SBT have been relatively low throughout the period, consistent with a high proportion of zero catches of SBT (Fig. S2), suggesting little or no targeting of SBT by Korean tuna longline vessels in these areas.

Figure 5 shows spatial patterns of SBT and ALB as proportions of the total reported catch south of 30°S by 5-year period. The SBT proportion was high in all periods, increasing further south, but declined in all areas after 2005. In the 2010s, particularly, there was little SBT taken in statistical area 8 north of about 37°S, whereas a high proportion of the catch in this area was ALB. This appears to reflect spatially and temporally differentiated targeting in area 8.

Sixty-five Korean tuna longline vessels have participated in statistical areas 8 and 9 since 1996 (Fig. S3), with over half of the fleet reporting their first participation before 2000. New vessels have arrived slowly but regularly. Vessel turnover was initially high, with over 20 vessels having stopped participating by 1998. In 2010, many vessels stopped participating because the SBT stock was at a critical stage, about 5% or less of the unfished spawning biomass level, and the quota was greatly reduced (CCSBT, 2009a; CCSBT, 2009b), and since then eight more have stopped but seven others have joined the fishery.

The total number of major cells (5° × 5° × month) fished has varied annually but declined considerably since the peak in 2009 (Fig. 6A). Over the same period, effort has become more concentrated with more operations per cell. This increasing concentration is also apparent at the minor cell (1° × 1° × month) level (Fig. 6B). The distribution of effort within major cells was more stable until recently, with similar numbers of minor cells per major cell on average, but in 2017-2018 increased to the highest level yet seen (Fig. 6C). Since 2008 the timing of effort in statistical areas 8 and 9 has changed, gradually moving earlier in the year, though with different timing peaks in each area.

Target change

Clustering

Applying Ward’s D hierarchical cluster analysis at the vessel-month level identified strong separation among two to three groups in statistical areas 8 and 9 (Fig. 7), so three clusters were chosen in each area to consider major targeting strategies. We preferred to use more clusters where there was uncertainty because unresolved target change can cause bias in indices.

In statistical area 8 (Figs. 8, 9 and 10), the species composition of cluster 2 was dominated by SBT, with small amounts of other species groups. Cluster 3 included similar amount of SBT, SHA and OTH, while cluster 1 included more ALB than SBT, with some OTH, YFT and BET. The SBT cluster 2 dominated the early part of the time series, with clusters 1 and 3 more apparent after 2005. Cluster 1 was represented during March to June, while clusters 2 and 3 occurred mostly in the second half of the year. Cluster 2 was fewer HBF than the average, while the hooks per set were similar for all clusters. Cluster 2 was well represented across most of the fished area, while cluster 3 occurred at middle latitudes from about 38°S–42°S, and cluster 1 occurred almost entirely in the far north of the area.

In statistical area 9 (Figs. 8–10), cluster 1 comprised almost entirely SBT, with small amounts of ALB and BET. Cluster 2 included significant SBT along with ALB, and some BET and YFT. Cluster 3 included similar amounts of SBT and OTH, with some SHA and ALB. Clusters 1, 2, and 3 were more strongly represented in the early, middle, and later parts of the time series, respectively. Clusters 1 and 3 occurred mostly in the period before August, while cluster 2 extended into October. The mean number of hooks was higher in cluster 3 and lower in cluster 2, and cluster 3 also had slightly more HBF. Cluster 2 dominated the northeast of area 9, while clusters 1 and 3 dominated the southeast and the southwest, respectively.

In summary, results show that effort in clusters 1 and 2 targeted ALB in area 8 and area 9, respectively. ALB targeting clusters operated further north than those targeting SBT. The main fishing periods for the ALB clusters were before June in area 8 and after June in area 9, whereas SBT targeting occurred after and before June, respectively.

CPUE standardization

All explanatory variables in the lognormal constant models and the lognormal components of the delta lognormal were statistically significant based on AIC (Table 1), with the year, location (latlong), vessel (vessid), and month effects the most important. The cluster effect was also important in statistical area 8, but less so in area 9. Several variables in the binomial component with weak support were retained in the model for consistency with the lognormal component, and because they had little effect on the results given the high rates of nonzero catches.

Nominal and standardized CPUE indices were developed for SBT in statistical areas 8 and 9, based on lognormal constant and delta lognormal models using data selection and cluster analysis (Fig. 11). The two methods to address targeting led to very similar standardized indices, with small differences from the nominal CPUE trend. Diagnostic frequency distributions and QQ-plots suggest that the data fitted the GLMs adequately (Figs. S5 and S6).

Differences between the targeting analysis methods were similar for both the lognormal constant and delta lognormal indices, though slightly larger for the lognormal constant method. The main differences between the methods occurred in the late 2000s for area 9 and in the 2013-2014 period for area 8, when indices were lower for the clustered data than the selected data. This may be because delta lognormal models are better than lognormal constant models at dealing with zero catches.

In addition, the area 9 delta lognormal indices differed from the lognormal constant indices. They were markedly lower before 2005 and were considerably higher in 2015 and 2018 but had similar trends to the nominal CPUEs in the recent years.

Hence the indices provided by the delta lognormal indices using the clustered data were chosen as the representative indices for SBT caught by the Korean tuna longline fishery. In summary, patterns in the indices differ somewhat between the statistical area 8 and 9 (Fig. 11B). Both sets of indices decreased until the mid-2000s, and subsequently increased, particularly in the last few years. However, lack of data prevents estimation for area 8 in the periods 2003-2007 and since 2017. Recent effort in area 8 is too low and concentrated to provide reliable estimates.

Influence plots for each covariate in the lognormal constant model using the clustered data are presented in Fig. 12, and those for the binomial and lognormal positive components in the delta lognormal model are presented in Figs. S7 and S8. Each subplot has three components, with the parameter estimates for each covariate level at the top, the effort by time interval and covariate level indicated by circle areas in the lower left component, and the cumulative influence of the covariate on the year effect on the right. Each subplot reports influence on a different scale, so the scales must be considered when comparing the relative importance of each covariate.

Vessel effects (Fig. 12A) were quite variable, with a few vessels having significantly lower or higher SBT catch rates. The influence changed through the time series, with the low number of vessels causing significant variability.

Spatial effects (Fig. 12B) showed significant variation in catch rates, with more variation in the statistical area 9 than area 8. In area 9 there was a trend towards fishing in areas with lower average catch rates. It would be useful to explore whether the areas of highest catch rate have moved through time. However, this would be difficult to determine from Korean data alone, since fishing activity is currently very concentrated spatially (Fig. 2). Given the behaviour of the species, areas of highest catch rate are also likely to move during the year, and we do not account for this in the model.

The effects of the number of hooks per set on catch rates (Fig. 12C) were difficult to interpret. In area 8 there were small differences by hook number across the range of data with most hooks, and minimal influence on year effects, apart from 2016 when effort was low and localized with only one vessel fishing. In area 9 there were relatively larger differences, and apparent influence on the year effects, with catchability averaging about 3% above the mean in 2010-16 (Fig. 12C). Sets with more than about 3,250 hooks tended to catch more SBT than sets with fewer hooks. In area 9 there were more sets with fewer hooks between 2004 and 2007, a period during which there were more zero SBT sets than at most other times (Fig. S2). These may reflect a mixture of targeting methods in area 9, with different fishing methods using different numbers of hooks.

Month effects (Fig. 12D) were strong in both areas 8 and 9, and relatively influential. In both areas, the highest catch rates were obtained in July and August, but slightly later in area 9. The seasonality of fishing effort changed through time, with the model suggesting that fishery timing increased mean catchability in the area 8 by over 10% higher than the average in 2015, and reduced it almost 5% below average in the area 9 in 2010-15.

Moon effects (Fig. 12E) showed that catch rates appeared to vary moderately with lunar illumination in area 8. In the area 9 delta lognormal models a similar effect was apparent in the lognormal positive component, but diminished overall by an inconclusive pattern in the binomial component (Fig. S8E).

The distribution of the cluster variable (Fig. 12F) changed considerably though time, as the behaviour of the fleet changed with the abundance of the target species. There were also relatively large differences in catch rate between the clusters, so this variable was quite influential on the indices, particularly in area 8.