This study investigated the temporal variability of YFT and BET catch amounts in the northern region (the Republic of Seychelles) and the southern region (the Republic of Mauritius + Tromelin Island + Reunion Island) of the SWTIO (Figure 1A). The annual catch amounts of the four EEZs were obtained from Sea Around Us (SAU, http://www.seaaroundus.org) for 1990–2014 (Boistol et al., 2011; Le Manach et al., 2015a; Le Manach et al., 2015b; Le Manach and Pauly, 2015). The catch data were based on officially reported data from the Food and Agriculture Organization of the United Nations (http://www.fao.org/fishery/statistics/en), but were improved by including unreported catches and major discards (Pauly and Zeller, 2015; Pauly and Zeller, 2016). However, the unreported catch amounts of YFT and BET are more than one order smaller than the reported catch amounts in the four EEZs. For example, the unreported YFT catch amount in the Republic of Seychelles is about 15 tonnes on average during 1990–2014.

The monthly sea surface height and upper 400 m temperature, salinity, and current velocity were obtained from the Simple Ocean Data Assimilation (SODA) 3.4.2, with a spatial resolution of 0.5° × 0.5° for 1988–2016 (Carton et al., 2018). The relationship between the tuna catches and the ocean environment was analyzed for the overlapping period, that is, 1990–2014, and the extra years, i.e., 1998–1999 and 2015–2016, were used to examine the lead-lag relationship with the climate indices. The thermocline depth and D20 were calculated from the SODA 3.4.2 data at each grid point, after conducting vertical linear interpolation of the temperature every 1 m. Monthly mean 10 m wind velocity and total precipitation were obtained for the same period (1988–2016) with a horizontal resolution of 0.25° × 0.25° from ERA5, which is the fifth generation atmospheric reanalysis from the European Centre for Medium Range Weather Forecasts. The ERA5 is based on 4D-Var data assimilation using global observation data (Hersbach et al., 2020). Monthly mean net primary production and chlorophyll concentration in the upper 400 m of the ocean were obtained for 1993–2016 with a resolution of 0.25° × 0.25°V from the Global Ocean Biogeochemistry Hindcast from Mercator-Ocean, provided by the Copernicus Marine Environment Monitoring System Service. This biogeochemical hindcast is based on the PISCES biogeochemical model, available on the Nucleus for European Modelling of the Ocean modelling platform (Perruche et al., 2019).

The CSEOF analysis (Kim et al., 1996; Kim and North, 1997) was used to extract the spatio-temporal variability of the environmental variables in the wider tropical Indian Ocean (30°N to 30°S and 40°E to 100°E, Figure 1A), including sea surface height, temperature, salinity, current velocity, 10 m wind velocity, total precipitation, net primary production, and chlorophyll concentration. The CSEOF technique has been applied by various studies to understand the relationship between a specific time series and the spatio-temporal variability of geophysical variables (Kim et al., 2015), including a study on the relationship between tuna fisheries and ocean environmental variability in the western equatorial Pacific (Kim et al., 2020).

In CSEOF, the environmental variables at each depth, A(r, t), were decomposed as cyclostationary loading vectors (CSLVs), B_n (r, t)VV and their corresponding principal component (PC) time series, C_n (t) as in Eq. (1):

where n, r, and t denote the mode number, 2-dimensional space, and time, respectively. Note that the CSEOF analysis was separately conducted on each environmental variable at each depth. For example, a total of 22 CSEOF decomposition was carried out for temperature, as the SODA provides 22 vertical levels in the upper 400 m. CSLVs illustrate the physical evolution of environmental variables with a nested period d, which was set to 12 months in this study. The nth mode of the CSLVs is expressed with periodicity d, as in Eq. (2):

The corresponding PC time series, C_n (t), illustrates the temporal variation of B_n (r, t) for 1988–2016, except for net primary production and chlorophyll concentration, which are available for 1993–2016.

A regression analysis was applied for each environmental variable at each depth, targeting the annual fishery time series, F(t) of the northern and southern regions separately, as shown in Eq. (3):

where m denotes the mode number, α_m is the regression coefficient, and is the annually averaged mth mode of the PC time series extracted from the environmental variables (predictor variables). The regression was conducted for the period of the target fishery data, 1990–2014 (1993–2014 for net primary production and chlorophyll), using the annual mean of C_m (t), because the target time series, F(t), are only available as the annual mean catch amounts. Note that the predictor PC time series, C_m (t), are uncorrelated to each other for each regression.

The regressed and reconstructed CSLVs, R(r,t), that illustrate the environmental variability related to the annual catch amounts, are obtained by using the regression coefficients, α_m as in Eq. (4):

where each R(r, t) consists of 12 monthly spatial patterns of each variable at each depth, same as B_m (r, t). Thus, the R(r, t) exhibits consistent temporal variability with the target time series F(t) with the error in Eq. (3). The number of PC time series used in the regression was set to ten in order to avoid a possible overfitting or underfitting issue. The R-squared values of the regression are shown in Supplementary Figure 1, and those of net primary production and chlorophyll concentration are shown in Supplementary Figure 2. If the R-squared value is relatively low for a specific variable at a specific depth, it indicates that the ocean environmental variability at the depth is not closely related to the target variability, i.e., the annual catch amount of YFT or BET in the northern or southern regions of the SWTIO. The annually averaged spatial patterns of each R(r, t)are presented in this study, focusing on the interannual variability.

The relationship between the ocean environmental variability and the IOD and ENSO was investigated using the same regression analysis, but using the monthly PC time series of the environmental variables, as the Indian Ocean Dipole Mode Index (DMI) and the Multivariate ENSO Index (MEI) are available as monthly means (Saji et al., 1999; Zhang et al., 2019). The regression analysis was also conducted by employing time lags between the ocean environmental variables and climate indices because the IOD and ENSO exhibit low-frequency variabilities that last longer than a year. A series of lagged regression analyses were conducted with time lags of ± 1 year (e.g., 1989–2013 and 1991–2015) and ±2 years (e.g., 1988–2012 and 1992–2016), targeting the DMI and MEI. The R-squared values of the lagged regression analysis are shown in Supplementary Figure 5. The R-squared values of the regression are shown at each depth in Supplementary Figure 3.

Figure 2 illustrates the zonally averaged mean and standard deviation of the temperature, salinity, and zonal and meridional current velocity in the upper 400 m of the SWTIO. The mean thermocline depth is shallower in the northern region (northern: ~66 m, southern: ~87 m), which is associated with the SCTR, and the peak of the SCTR (shallowest thermocline depth) is observed at approximately 8°S (Figure 2A). The mean D20 is shallower (~91 m) and closer to the thermocline depth in the northern region, but the mean D20 is substantially deeper (~163 m) in the southern region. In the northern region, the near-surface (< D20) is more saline, but the deeper layer (> D20) is more saline in the southern region (Figure 2B). The fresh near-surface water in the southern region is from the Indonesian Throughflow in the equatorial eastern Indian Ocean (Gordon et al., 1997; Makarim et al., 2019). The mean zonal currents show meridionally opposing directions, reflecting the southeastward-flowing South Equatorial Counter Current (SECC) to the north of the SCTR (north of approximately 8°S) and the southwestward-flowing South Equatorial Current (SEC) to the south of the SCTR (Figures 2C, D).

The temporal variations in the environmental variables also exhibit meridional differences on the interannual time scale (see the bottom of Figure 2). Due to the strong subsurface upwelling variability in the SCTR, the standard deviation of the upper ocean (< 400 m) temperature is larger near the thermocline depth and D20 in the northern region. In both regions, salinity exhibits significant variability from the surface to the thermocline depth, which is deeper in the southern region. The standard deviations of the zonal and meridional currents are larger in the northern region, where the SECC flows southeastward. The different ocean environmental conditions in the northern and southern regions shown in Figure 2 may contribute to the YFT and BET catch variability in the SWTIO.

Figure 3 displays the annual catch amount of the YFT and the BET in the SWTIO for 1990–2014. The YFT and BET catch amounts throughout the SWTIO exhibit similar variability, except for a few years in the early 2000s (Figure 3A). This similarity between the YFT and BET catch amounts is also observed separately in the northern and southern regions (Figures 3B, C). In the northern region, both YFT and BET catches were larger in the 1990s and the early 2000s, compared to the late 2000s and early 2010s. In the southern region, however, the catches were smaller in the 1990s than in later years. The difference between the YFT and BET catch variability seems not to be significant as the difference of the catch variability between the northern and southern regions.

Figure 4 illustrates the regressed anomalies of the ocean conditions during the increase of the YFT and the BET in the northern and southern regions of the SWTIO. The R-squared values of the regression for each variable and each depth are presented in Supplementary Figure 1. The temperature conditions, during the increase of the YFT and the BET catches in the northern region, represent strong positive anomalies near the mean thermocline depth and D20 in the northern region (top two rows in Figure 4A). They also exhibit the negative anomalies below the thermocline depth and D20 in the northern SCTR. The temperature anomalies indicate the deeper isotherms in the southern SCTR and the shallower isotherms in the northern SCTR, thus the flattening of the isotherms in the SCTR, during the increase of the catch in the northern region. The flattening of the isotherms suggests the suppression of the subsurface upwelling in the SCTR. On the contrary, the temperature conditions, during the increase of the YFT and the BET catches in the southern region, represent relatively weak positive anomalies in the southern region (bottom two rows in Figure 4A).

The salinity conditions, during the increase of the YFT and the BET catches in the northern region, exhibit similar anomalies with the near-surface (< D20) freshening in the northern region (top two rows in Figure 4B). Those freshening in the northern region is consistent with an increase of precipitation in this region (top row in Supplementary Figure 4). On the other hand, the positive salinity anomalies and the decrease of precipitation are detected in the southern region during the increase of the YFT and the BET catches in the southern region (bottom two rows in Figure 4B and middle row in Supplementary Figure 4). The negative anomalies in the north and the positive anomalies in the south contribute to the similar near-surface salinity ranges of 34.8–35.2 psu, which is within the favorable salinity range for an increase of the YFT and the BET catches reported in the Indian Ocean (Song et al., 2008; Song et al., 2009).

The current velocity conditions in the northern region during the increase of the YFT and the BET catches in the northern region exhibit larger anomalies of zonal current velocity compared to those of meridional current velocity (close to an order of magnitude of difference) (top two rows in Figures 4C, D). They exhibit westward anomalies at the latitude near the SCTR upwelling peak and the eastward anomalies to the north of the peak. It indicates the equatorward shift of the SEC in comparison with the mean current velocity in Figures 2C, D. The equatorward shift is consistent with the easterly wind anomalies near the equator during the increase of YFT and BET catches in the northern region (top row in Supplementary Figure 4). The regressed anomalies of current velocity targeting the YFT and the BET catches in the southern region are relatively smaller (bottom two rows in Figures 4C, D), similar to the results of the temperature regression.

Figure 5 illustrates the ocean environmental conditions in the SWTIO during the positive phase of the IOD and the ENSO. The R-squared values of the regressions are displayed in Supplementary Figure 3. The strong positive temperature anomalies are observed near the thermocline depth and D20 in the northern region during the positive phase of both the IOD and the ENSO, indicating the suppression of SCTR upwelling (Figure 5A). The current velocity conditions represent an equatorward shift in the SEC during the positive phase of both the IOD and the ENSO (Figures 5C, D).

However, the salinity conditions are different between the positive phase of the IOD and the ENSO (Figure 5B). During the positive phase of the IOD, the salinity anomalies exhibit different signs between the northern and southern regions (top of Figure 5B). The dipole-like pattern of salinity anomaly was previously reported as a Salinity-IOD Mode that originated from the eastern tropical Indian Ocean (Zhang et al., 2016). On the other hand, during the positive phase of the ENSO, the salinity condition reflects an overall near-surface freshening in both northern and southern regions (bottom of Figure 5B). The regression of the precipitation generally exhibits consistent patterns with the regressed salinity anomalies (Supplementary Figure 4).