The weighted habitat suitability index (HSI) was calculated using environmental factors – the SST, bottom temperature (BT), sea surface height (SSH), chlorophyll-a (Chl-a), and salinity. The HSI model coupled with a climatic and environmental variability analysis can provide insights into climate-driven habitat suitability variation under different anomalous events for swimming crab species. To elucidate the mechanisms through which climate change affects crab habitat suitability, oceanographic characteristics and suitable habitat areas for the three swimming crabs in the TS during various ENSO events were compared (Figure 1).

Fishery data on traps from 2013 to 2019 were collected from the voyage data records and logbooks of 218 Taiwanese crab vessels (ranges from 61 to 118 vessels across years with 5–200 metric tons) operating in the TS (Supplementary Table 1 and Supplementary Figure 1). The fishing data comprised daily fishing positions for 0.1° spatial grids, including latitude and longitude, fishing date, soaking time, total catch (kg), and species. C. feriatus, P. pelagicus, and P. sanguinolentus were selected for analysis and other species recorded were excluded. Swimming crabs were collected using circular crab traps with frozen mackerel as bait placed in the center of the trap.

Monthly CRs were calculated as the cumulative weight of each crab species caught from all fishing vessels for a month divided by the cumulative soaking time of all set lines in that month to adjust for fishing effort. Monthly CRs were calculated for individual 0.1° grid squares across the study region by using the following expression:

where i represents the latitude and j represents the longitude of each 0.1° spatial grid square.

We constructed the distribution maps of CRs for the three swimming crab species using Interface Descriptive Language version 7.0.

To obtain tropical climate signals, we considered the Oceanic Niño Index (ONI), an ENSO-related index and examined the relationship between the mean CRs of C. feriatus, P. pelagicus, and P. sanguinolentus. The ONI is currently used by National Oceanic and Atmospheric Administration to provide forecasts of El Niño and La Niña activity and is defined as the 3-month running mean of SST anomalies in the Niño 3.4 region of the equatorial Pacific (5°N–5°S and 120°–170°W). An El Niño and La Niña event is defined to occur when the ONI index higher than +0.5°C and lower than -0.5°C for 5 consecutive months, respectively. An ONI in the range ±0.5°C is defined as a normal year. The definition and the ONI data were sourced from the NOAA climate prediction center¹.

Marine environmental data – daily SST, BT, SSH, and sea surface salinity (SSS) – were extracted from the Copernicus Marine Environment Monitoring Service² from 2013 to 2019 at a horizontal resolution of 1/12°. The SST, SSH, and SSS were obtained across the upper 0.5 m of surface water and potential temperature at the seafloor of each grid for the BT.

Monthly Chl-a data were downloaded from the NASA Aqua satellite³, which has a sensor that detects the concentration of Chl-a in the world’s oceans. The level 3 data map image of Chl-a data has a monthly temporal resolution and a spatial resolution of 4.6 km (at the equator).

The daily environmental data were then calculated mean monthly on a spatial grid of resolution 0.1° to fit the fishery data. Intercorrelations among environmental data were tested before application of statistical models to determine the correlation coefficients among the variables. The test for intercorrelations among all environmental variables obtained low correlation coefficients (r² < 0.25).

The relationship of each swimming crab CRs with the environmental variables in ENSO events on the basis of crab trap fishery data locations from 2013 to 2019 in the TS were analyzed using a one-way ANOVA. Pairwise post hoc comparisons were conducted using Tukey’s tests. For the one-way ANOVA, we tested and found that required assumptions were met with regard to the homogeneity of variance, normal distribution of errors, and independence of errors.

The habitat suitability modeling approach helps efficiently obtain the distribution pattern of potential habitats for marine species, including crustaceans (Chen et al., 2019). The HSI is a numerical index between 0 and 1, where 0 indicates an unsuitable habitat and 1 represents an optimal habitat. Suitable habitat is closely related to the area that had the highest catches (Chang et al., 2019). By considering CRs > 60% as the high catches, the relationships between the CRs for three swimming crabs and each environmental variable were presented as curves of the suitability index (SI, 0 to 1). The formula for the SI was as follows:

where catch rates_ij denotes the cumulative fishing efforts in the ith interval of the range of environmental variable j, and max(catch rates_ij) indicates the maximum cumulative fishing efforts. The relationships between the SI and environmental variables for each swimming crab species were determined as follows (Yu et al., 2018):

where a and b are the model parameters to be estimated, and i denotes the environmental variables considered in this study.

The fitted SI models were combined into an integrated HSI model by using the arithmetic mean model; this approach is the most commonly used and yields optimal performance (Silva et al., 2016; Yu et al., 2018). Studies have reported that setting weights for different environmental variables can significantly influence HSI modeling.

Partial least squares regression (PLSR) is a bilinear calibration method that compresses data by reducing many measured collinear spectral variables to a few non-correlated principal components. PLSR is a valuable tool for exploring datasets with numerous predictor variables (Carrascal et al., 2009). PLSR considers the latent structure of the predictor variables (CRs > 60%) and response variables (five environmental parameters) and selects the model with predictor variables that best describe the variation in the response variables (Wold et al., 2001). PLSR variable influence on projection (VIP) scores were generated to assess the relative importance of the predictors in the overall model, and the relative weighting was based proportionally on the VIP score.

Therefore, in this study, five weights based on VIP scores were used for each swimming crab species in the HSI model development. The weighted HSI arithmetic mean model can be described as follows:

where SI_i is the suitability index or value associated with environmental variable i, n is the number of environmental variables considered in this study, and w_i is the weight based on the VIP scores of each environmental variable for each swimming crab species. The spatial distribution of the weighted HSI obtained from the models in the major fishing seasons was mapped, and the influence of ENSO events to habitat suitability of C. feriatus, P. pelagicus, and P. sanguinolentus was investigated.

The monthly mean catch and CRs in the TS for each species exhibited similar trends (Figures 2A–C). The highest monthly mean catch and CRs for the three swimming crabs were observed from September to November (autumn) and then gradually decreased until June to July. Because the major fishing season with the highest catch and CRs for C. feriatus, P. pelagicus, and P. sanguinolentus was the autumn, this season was selected as the study period to determine the effects of climate events.

The CR distribution maps for C. feriatus revealed that higher CRs were distributed in the northern part of the TS from 25 to 26°N and 120 and 122°E during autumn (Figure 3A). For P. pelagicus, higher CRs were distributed further to the south of the TS between 22 and 24°N and extended to the northeast at 25–26°N (Figure 3B). P. sanguinolentus was widely distributed across the TS, with higher CRs extending to northeastern areas during autumn (Figure 3C).

Based on the definition of ENSO events, La Niña events occurred during the autumns of 2016 and 2017, and El Niño events occurred during the autumns of 2014, 2015, and 2018 (Figure 4). Compared with the mean CRs of ENSO events, the CRs for C. feriatus and P. pelagicus were higher (>7.0 and >8.0 kg/h) in La Niña events and were more than 40.0% higher compared with the CRs in normal and El Niño events during autumn (Figure 4). For P. sanguinolentus, the effect of ENSO events was not clear; the CRs were higher in both La Niña and El Niño events (>8.0 kg/h) compared with those in normal years.

The mean CRs were further calculated for the northern TS, middle TS, and southern TS (Figure 3) to understand the effect of ENSO events on the distribution of three swimming crabs during autumn (Figure 5). The CRs for C. feriatus were higher in the northern TS than in the middle or southern TS, especially in 2017, a La Niña year (Figure 5A). The lowest CRs for C. feriatus in the northern TS occurred in the El Niño events of 2015, and the CRs were lower than 1.7 kg/h in the southern TS throughout the study period. For P. pelagicus, the mean CRs were higher in the southern TS and lower in the middle and northern TS (Figure 5B). In the southern TS, the highest CR for P. pelagicus was recorded in the La Niña year 2016, whereas the lowest CRs occurred in the normal year 2019. CRs < 1.0 kg/h were recorded for P. pelagicus in middle and northern TS for each climate event from 2013 to 2019. The CRs for P. sanguinolentus were higher in the northern and southern TS throughout the study period than in the middle TS, especially in 2017, a La Niña year (Figure 5C). The lowest CRs for P. sanguinolentus (<3 kg/h) occurred in the middle TS during the normal year 2019.

Obtained through PLSR, the relative weighting of important environmental variables of the CRs for C. feriatus, P. pelagicus, and P. sanguinolentus based on VIP scores are presented in Table 1. The BT explained the highest VIP score for C. feriatus, followed by the SSH, SSS, and Chl-a, whereas the SST explained the lowest (Table 1). Similarly, for P. sanguinolentus, the BT explained the highest VIP score, followed by the Chl-a, SSH, SST, and SSS. By contrast, for P. pelagicus, SSS explained the highest VIP score, followed by the SST, BT, and SSH, with Chl-a explaining the lowest score.

The box-and-whisker plots illustrate the relationship of each swimming crab species with the environmental variables during autumn in ENSO events on the basis of crab trap fishery data locations from 2013 to 2019 in the TS (Figure 6). Statistically, ANOVA showed significant effects of environmental variables during autumn in ENSO events on CRs of all swimming crab species (Figure 6). The results revealed that high CRs for C. feriatus in La Niña co-occurred with higher SST, BT, and SSH and lower Chl-a and SSS (Figures 6A–F). However, lower CRs for C. feriatus in normal and El Niño years revealed no significant differences and co-occurred with lower SST, BT, and SSH and higher Chl-a.

For P. pelagicus, the CRs were higher in La Niña years and co-occurred with higher SST, BT, SSH, SSS, and Chl-a (Figures 6G–L). Similar to C. feriatus, no significant differences were noted in lower CRs for P. pelagicus in the normal and El Niño years, and they co-occurred with lower SST, SSH, SSS, and Chl-a. Moreover, for P. sanguinolentus, higher CRs were recorded in La Niña and El Niño years and co-occurred with higher BT and SSH, whereas the SST and Chl-a were higher in La Niña years and the SSS was not significantly different between various ENSO events (Figures 6M–R).

The fitted SI model was established for environmental variables in autumn, and all the models were significant (p < 0.01) with low root-mean-square errors (RMSEs) and high correlation coefficients (Table 2). Different optimal ranges for each environmental variable were identified for each swimming crab species when the SI curve was used (Figure 7). The optimal range for SSTs with higher SI (>0.6) for C. feriatus and P. sanguinolentus varied from approximately 25.0 to 29.0°C, whereas that for P. pelagicus varied from 27.0 to 30.0°C (Figure 7A). For BTs, the optimal ranges were approximately from 19.0 to 26.0°C for C. feriatus and P. sanguinolentus and from 24.0 to 27.0°C for P. pelagicus (Figure 7B). The optimal range for SSHs for C. feriatus and P. sanguinolentus was approximately from 0.52 to 0.64 m, whereas that for P. pelagicus was from 0.64 to 0.76 m (Figure 7C). Moreover, the optimal range for SSSs was 32.8–33.7 PSU for C. feriatus and P. sanguinolentus and 33.1–34.0 PSU for P. pelagicus (Figure 7D). However, the suitable Chl-a range was higher for C. feriatus than for P. pelagicus and P. sanguinolentus (0.7–1.4 mg/m³; Figure 7E).

The study of HSI maps revealed that suitable habitat areas (HSI > 0.6) of C. feriatus in normal (2013) and El Niño (2014) years were smaller in the northern TS (Figures 8A,B). By contrast, the suitable habitat areas of C. feriatus were higher in the northern TS and extended to northern areas (north to 26°N) in La Niña year 2017 (Figure 8C), with higher CRs and a higher HSI, indicating that fishing locations moved further north in the TS.

For P. pelagicus, the selected habitat maps show suitable areas in the southern TS (south to 24°N) in each ENSO event (Figures 8D–F). However, the suitable habitat area was smaller during autumn in the normal (2013) and El Niño (2014) years than in the La Niña years. Moreover, a smaller HSI indicated that fishing locations moved to the northern TS (25–24°N) during the normal and El Niño years (Figures 8D,E).

The HSI maps reveal that the suitable habitat areas of P. sanguinolentus during autumn in a normal year (2013) were lower in the northern TS and further to the south of the TS (Figure 8G), whereas the suitable habitat area was higher in the northern TS during autumn in the La Niña years (2017) (Figure 8I). By contrast, the suitable habitat of P. sanguinolentus was widely distributed across the TS during autumn in El Niño years (Figure 8H).

A further comparison of the interannual variability of CRs with HSI in ENSO events (Figure 9) revealed that higher CRs for C. feriatus and P. pelagicus were distributed in the higher HSI areas between 25 and 26°N and between 22.5 and 23.5°N, respectively, in the La Niña years (Figures 9A,B). However, the higher CRs areas were distributed more widely further north of 26.5°N for C. feriatus and 23.5°N for P. pelagicus during normal and El Niño years. The high CRs for P. sanguinolentus were widely distributed in the TS in each climate event (Figure 9C).

This study demonstrated that autumn is a crucial season for obtaining large catches and high CRs for C. feriatus, P. pelagicus, and P. sanguinolentus in the TS. Furthermore, high CRs for P. sanguinolentus were widely distributed in the TS, whereas those for P. pelagicus and C. feriatus were found in two main areas in the southern and northern TS during autumn. The optimal season for crab fishing is dependent on various factors such as population size, food availability, relevant environmental factors, and the type of gear used (Olsen et al., 2019; Spencer et al., 2019; Naimullah et al., 2020). Variation in CRs is most likely due to environmental factors affecting recruitment success over time and the geographic distribution of swimming crabs in the TS.

The decrease in catches and CRs for C. feriatus, P. pelagicus, and P. sanguinolentus from March to July may have been related to the spawning seasons and recruitment process. Brachyuran crabs inhabiting tropical waters usually breed throughout the year, whereas those found in temperate waters breed only in certain months (Warner, 1977). The spawning season of the swimming crab in the East China Sea was recorded from March to July (Song et al., 1989; Wu et al., 2007). Moreover, because of low water temperatures during winter and spring in the East China Sea, berried females do not hatch until late April (Song et al., 1989; Xie et al., 2002). The berried swimming crab mostly remains buried in the sand and carries its eggs until the larvae hatch as zoeae, which drift about in the water current (Rasheed and Mustaquim, 2010). During spawning seasons, female swimming crabs have a period of low catchability when they cease feeding and become inactive in preparation for the molt and are not easily catchable with baited traps or drop nets (Soundarapandian et al., 2013). Therefore, we suggest that autumn is the preferred season for fishing of these three swimming crab species in the TS.

The higher CRs for P. pelagicus and C. feriatus in the southern and northern TS during autumn may have been due to the habitat and environment being conducive to spawning, recruitment, and growth. After the hatched zoea transform into juvenile crabs, most of those in coastal habitats settle there or migrate to estuaries as conditions in the estuaries become favorable. Moreover, some studies have shown that swimming crabs that inhabit marine embayments do not leave these marine environments to spawn (Potter and De Lestang, 2000; Polity et al., 2011). The wide distribution of P. sanguinolentus in the TS may be primarily influenced by sea currents. Most P. sanguinolentus larvae hatch as zoea, which drift in the water current before inhabiting benthic habitats and growing into juvenile crabs (Rasheed and Mustaquim, 2010). Throughout the year, the strong Kuroshio Branch Current and South China Sea Current in the TS during autumn may affect the movement of P. sanguinolentus in the northern area. P. sanguinolentus typically inhabit the sandy and muddy bottom in shallow areas (Sumpton et al., 1989; Carpenter et al., 1997), whereas adult and berried females often migrate to deeper waters for spawning (Yang et al., 2014), which results in a widespread distribution of this species.

Understanding the relationships between environmental variables and fundamental processes is crucial to explaining variation in CRs. Although many environmental variables affected the abundance of the three investigated swimming crabs in the coastal waters of Taiwan (Naimullah et al., 2020), each factor contributed differently to the stock dynamics. BT had a positive impact (>21.7°C) on the habitat suitability for C. feriatus and P. sanguinolentus in autumn, implying that rising water temperatures are favorable to the formation of potentially high-quality habitats. Temperature is a crucial factor influencing the physiological and ecological properties of brachyuran crabs (Azra et al., 2020). The variation in water temperature influences molt timing and frequency, growth rates, size at maturity, spawning period, larval development, geographic range, and survival (Defeo and Cardoso, 2002; Ikhwanuddin et al., 2012; Green et al., 2014; Ryer et al., 2016; Waiho et al., 2016; Chandrapavan et al., 2019).

Water temperature is known to modify reproduction and maturation in crustaceans, particularly crabs (Olson et al., 2018). Cooler temperatures generally retard growth and delay maturity, causing animals, including crabs, to begin maturation when they are larger (Leffler, 1972; Azra et al., 2020). By contrast, rearing crabs at higher water temperatures may benefit their growth, molting, and time to maturity (Kuhn and Darnell, 2019). Reducing the size of mature female crabs may negatively affect reproduction because body size is a factor limiting fertility (Green et al., 2014). Environmental variations may retard brachyuran crab growth (Azra et al., 2020). High CRs for C. feriatus and P. sanguinolentus in autumn were also found when the BT was higher during La Niña years, which affected the catch and life cycle processes of these species. Moreover, the La Niña event in Northern Territory Australia showed that intense rainfall and warm waters may explain between 30 and 40% of the Scylla serrata catch variability (Meynecke et al., 2012a).

The SSS, in the range 33.1–34.0 PSU, had a positive impact on habitat quality for P. pelagicus and was higher in La Niña years in the fishing positions. Batoy et al. (1987) also discovered that salinity was one of the most critical environmental factors affecting the reproductive cycle of P. pelagicus species, with higher salinity more favorable for breeding. Furthermore, the optimal range of the SSS was 32.8–33.7 PSU for C. feriatus and P. sanguinolentus in the TS. In Haizhou Bay, salinity significantly influenced the distribution of C. bimaculata, indicating that the preferred salinity range was 29–31 PSU (Luan et al., 2018). By contrast, P. sanguinolentus females mainly inhabit deeper waters and prefer higher salinity than male crabs (Campbell, 1986). Many crustacean fisheries are located in coastal and estuarine areas where the salinity fluctuates due to rainfall and freshwater in-flow (Meynecke et al., 2012a,b). Crabs are generally restricted by a lower salinity threshold; therefore, both adults and larvae avoid low salinity by modifying their behavior such as migration (Gibson and Najjar, 2000; Charmantier et al., 2001).

The Chl-a concentration explained one of the vital environmental factors for P. sanguinolentus. Signa et al. (2008) determined that the density of the swimming crab Polybius henslowii is strongly related to the Chl-a concentration, suggesting higher accumulation in locations with higher production. Moreover, the highest CRs for P. sanguinolentus in the TS were mostly in waters with Chl-a concentration >0.5 mg/m³ (Naimullah et al., 2020). The essential fishing ground was mainly related to the upwelling area consisting of a high Chl-a concentration and low SST that carried large amounts of nutrients from the bottom layer to the ocean’s surface layer (Hu and Wang, 2016) and is a crucial factor affecting the diet of P. sanguinolentus in TS. The SSH, a proxy for the current and temperature balance, was one of the critical factors affecting the swimming crab density. Liu et al. (2001) reported that changes in thermocline depth were related to the SSH. In addition, an increase in the SSH would increase the temperature, thus affecting the crab habitat. A greater SSH (>0.60 m) in La Niña years yielded a positive impact on the CRs for the three swimming crabs in the TS. Srisurichan et al. (2005) reported a similar outcome for spiny lobsters; they reported a positive response of CRs for legal-sized spiny lobsters to the SSH. They noted that swell might produce conditions favoring better protection from predators by increasing bottom turbidity and food availability.

A modeling approach for evaluating the quality of crustacean species’ habitats is attracting considerable attention and has been increasingly applied in fisheries research (Luan et al., 2018; Chen et al., 2019; Naimullah et al., 2020). The El Niño and La Niña events in the Pacific Ocean represent considerably complex climate variability (Wang et al., 2017). The present study found that ENSO events have substantial impacts on the environmental conditions of fishing grounds for C. feriatus, P. pelagicus, and P. sanguinolentus and change the suitability of habitats in the TS. Higher CRs were observed in La Niña years than in normal and El Niño years for C. feriatus and P. pelagicus, whereas no significant effect of ENSO events was noted on the CRs and distributions of P. sanguinolentus. Furthermore, higher CRs for C. feriatus and P. pelagicus were concentrated in the higher-HSI areas in the La Niña years but distributed more widely in areas with lower HSI during normal and El Niño years. These findings indicate that fishing vessels require more time and oil to search the fishing grounds for C. feriatus and P. pelagicus during normal and La Niña years.

Climate-induced habitat suitability changes may strongly influence fisheries (Kuczynski et al., 2018). For example, in the southern TS, pelagic species resources are susceptible to climate change, whereas benthic species are mainly insensitive to climatic factors (Hsiao et al., 2021). This variability may explain changes in CRs in suitable and optimal habitat areas caused by ENSO events in the TS. In each season, the Kuroshio Branch proportions are the highest during El Niño periods, and this leads to higher salinity than that during La Niña events in the TS (Huang et al., 2015). Concentrations of most nutrients, hydrogen ions, and Chl-a are higher in the autumn of a La Niña year than in the autumn of a normal or El Niño year in the southern TS (Huang et al., 2015). However, the SST was found to be lower in an El Niño summer than in a La Niña summer along the central TS northeastward and up to 24°N (Kuo and Ho, 2004). Thus, we suggest that higher nutrient levels and sea temperature lead to greater SSH in autumn and more extensive optimal habitat areas, higher CRs, and greater growth of C. feriatus and P. pelagicus in the La Niña year in the TS. However, the CRs for P. sanguinolentus are not influenced by ENSO events due to its wide distribution and optimal habitat areas and these species being able to strongly adapt to changing environments.

Numerous environmental factors affect crabs throughout their complex life cycle (Green et al., 2014), and lagged effects on swimming crabs have been reported (Meynecke et al., 2012a; Chandrapavan et al., 2019). Meynecke et al. (2012a) indicated that the best environmental predictor of S. serrata catch variability was the Southern Oscillation Index with a 6–9-month lag; negative impacts on catches were observed after heavy and prolonged flooding with a 1–2-year lag effect. However, in Western Australia, trap CRs and the monthly commercial catch from both trap and trawl fleets declined to historically low levels almost 6 months after extreme marine heatwave phenomena in the 2011 La Niña (Chandrapavan et al., 2019). In addition, swimming crabs mostly have a maximum lifespan of approximately 3 years, and the abundance of crabs is highly dependent on the success of recruitment and the survival of a small number of year classes (Sukumaran and Neelakantan, 1997). High fishing pressure on crab stocks, unsuitable habitats, and climate change year can reduce a year’s swimming crab catch (Meynecke et al., 2012a,b; Chandrapavan et al., 2019). Although our study only focused on the season in which swimming crab production is highest, we expect that negative or positive climate change events in previous seasons would influence the swimming crab CRs in the following season, which should be further studied.