Assessing Summertime Primary Production Required in Changed Marine Environments in Upwelling Ecosystems Around the Taiwan Bank

: The Taiwan Bank (TB) is located in the southern Taiwan Strait, where the marine environments are affected by South China Sea Warm Current and Kuroshio Branch Current in summer. The bottom water ﬂows upward along the edge of the continental shelf, forming an upwelling region that is an essential high-productivity ﬁshing ground. Using trophic dynamic theory, ﬁshery resources can be converted into primary production required (PPR) by primary production, which indicates the environmental tolerance of marine ecosystems. This study calculated the PPR of benthic and pelagic species, sea surface temperature (SST), upwelling size, and net primary production (NPP) to analyze ﬁshery resource structure and the spatial distribution of PPR in upwelling, non-upwelling, and thermal front (frontal) areas of the TB in summer. Pelagic species, predominated by those in the Scombridae, Carangidae families and Trachurus japonicus , accounted for 77% of PPR (67% of the total catch). The benthic species were dominated by Mene maculata and members of the Loliginidae family. The upwelling intensity was the strongest in June and weakest in August. Generalized additive models revealed that the benthic species PPR in frontal habitats had the highest deviance explained (28.5%). Moreover, frontal habitats were inﬂuenced by NPP, which was also the main factor affecting the PPR of benthic species in all three habitats. Pelagic species were affected by high NPP, as well as low SST and negative values of the multivariate El Niño–Southern Oscillation (ENSO) index in upwelling habitats (16.9%) and non-upwelling habitats (11.5%). The composition of pelagic species varied by habitat; this variation can be ascribed to impacts from the ENSO. No signiﬁcant differences were noted in benthic species composition. Overall, pelagic species resources are susceptible to climate change, whereas benthic species are mostly insensitive to climatic factors and are more affected by NPP. Author Contributions: Conceptualization, K.-W.L.; methodology, K.-W.L. and P.-Y.H.; software, T.S.; formal analysis, K.-W.L. and M.-A.L.; resources, C.-H.L.; data curation, P.-Y.H., M.-A.L. and C.-H.L.; writing—original draft preparation, P.-Y.H.; writing—review and editing, P.-Y.H. and K.-W.L.; supervision, K.-W.L.


Introduction
Located in the southern part of the Taiwan Strait (TS), the Taiwan Bank (TB) is characterized by sand dunes and shallow water (approximate depth 10-30 m). Three primary currents affect the waters around the TB: the Kuroshio Branch Current (KBC), South China Sea Warm Current (SCSC), and China Coastal Current (CCC). These currents are defined by the seasonal monsoons [1][2][3]. In summer, the northeastward surface currents of the SCSC and KBC bringing the warm, saline, and oligotrophic water from the SCS to the ECS, which are predominantly driven by wind, and the bottom currents flow upward from the continental slope, forming upwellings with vertical mixing of different layers' water ( Figure 1) [1][2][3][4][5]. The upwelling and thermal front regions near the southern edge of the TB span approximately 2500-3000 km 2 and are clearly observable in summer [2,4,5]. Thermal fronts with characteristic shapes appear at fixed locations around the upwelling, where the sea surface temperatures (SSTs) are lower and primary production (PP) is higher than those in the surrounding waters [5,6]. Here, four different water masses (the diluted water, the Kuroshio water, the upwelling water, and the warm water from the South China Sea (SCS)) contribute to the formation of fronts [5]. The thermal fronts and orientations around TB upwelling correspond to cold fronts with characteristic shapes in the north and south of TB upwelling [6]. Furthermore, small and mesoscale eddies, filaments caused by the monsoons and KBC, which all have an effect on the larva species due to the dispersion and convergence of nutrients [7,8]. Through foraging behavior, the aggregative response of biological interactions links primary productivity and top predators. Finally, the distribution of fishery resources in the different fisheries is affected [9][10][11][12]. High PP contributes to abundant fishery resources. The upwelling of nutrient-rich bottom water into the euphotic zone promotes the growth and propagation of phytoplankton, sustaining zooplankton communities and generally supporting fishery production. Therefore, upwelling strength directly affects the stability of marine ecosystems and general oceanic productivity [4,12,13]. Oceanic fronts (e.g., thermal fronts), sharp transition regions between water masses with disparate hydrographic, biological, or chemical properties, can limit the distribution of marine organisms, thereby serving as biogeographical boundaries [14][15][16]. The hydrographic characteristics of thermal fronts provide environmental conditions conducive to plankton blooms, fish aggregations, and consumer congregation [17,18]. Upwellings and thermal fronts have supported the maintenance of critical commercial fisheries near the TB since the 1950s; such fisheries typically employ light fishery, pole and line, and longline fishing, as well as trawling and gill nets [19]. The total catch weight and value of the offshore and coastal fisheries in the TS peaked in the 1980s, accounting for 44% of Taiwanese fishery production. By 2018, it had decreased to 16% ( Figure 2). This dramatic decline can be attributed to high demand, unregulated fishing practices, and overexploitation [20,21]. In their evaluation of the stock status of 16 species under fishing pressure in the waters surrounding Taiwan, Ju et al. [21] discovered only one species that is not overfished. Furthermore, climate-driven environmental variabilities such as the El Niño-Southern Oscillation (ENSO) are pivotal factors affecting the intensity of the upwelling regions in the TB. The East Asian summer monsoon is always weaker at the peak of El Niño but strengthens when El Niño starts to decay [22][23][24]. This suggests that stronger south westerlies occur in the summer following an El Niño event and in turn more strongly promote coastal upwelling. In addition to warming the SCS currents northward, these winds produce higher SST gradients in the TS [12,25]. Changes in marine environments have also been demonstrated to affect fishery resources on different spatiotemporal scales, as well as ecological factors such as ecotourism and fisheries landings [26][27][28]. The marine environmental factors not only directly affect the habitat environment of biological preference but are also impactful on the competition and predation relationship of species in the ecosystem [27]. The bottom-up and top-down control could limit all the taxa in the ecosystem with interaction between species and the energy flow [29]. PP, the basis for food chains, promotes the abundance of fishery resources [30][31][32]. Characterized by energy transfer between trophic levels, PP can be converted into percentages relative to the primary production required (PPR) to provide an estimate for the sustainability of fishery production [33]. In satellite telemetry studies, Chassot et al. [34] and Lam and Pauly [35] have found large marine ecosystems (LMEs) in tropical, temperate, and upwelling systems to have varying levels of fishery sustainability, trophic transfer efficiency, and species composition-all of which affect PPR. Watson et al. [36] investigated the historical operation of fishery fleets in major fishing grounds near six continents to gain an understanding of specific ecosystems and national fisheries, comparing their potential fishery sustainability. Knight and Jiang [37] used a fishing pressure index that compares NPP and PPR estimates for a specific catch to assess the sustainability of local fisheries in New Zealand and thereby provide a reference for subsequent fishery management. PPR can also be applied to large-scale fisheries, such as those in LMEs. Taken together, these studies indicate that PPR can reflect changes in fish catch that are related to energy flow; therefore, the interspecies transfer efficiency may accurately reflect the effect of ecosystems on fisheries.
Along with upwellings, the complexity of hydrological cycling and relevant changes that occur in summer are responsible for the high nutrient levels and high PP in the TB, a diverse fishing environment characterized by richness of species and habitats, which include upwelling, frontal, and non-upwelling habitats. However, studies on the TB have devoted little attention to differences in habitat characteristics when investigating the impacts of climate variability on PP and PPR in this region. Therefore, the present study collected data on fishery activity in the TB area and investigated the fishery resource structure and spatial distribution of PPR in upwelling, non-upwelling, and thermal front (frontal) habitats in summer. In addition, we explored the influence of marine climatic changes on the fishery resource structure and PPR of different habitat characteristics. The results are a valuable reference for ecosystem-based fishery management.

Fishery Data Collection
High spatial-resolution data on 11 types of fishery activity were collected from the voyage data recorders and logbooks of fishing vessels operating in the waters of the TB (118 • E-120 • E, 22 • N-24 • N) from 2011 to 2016 ( Table 1). The fishery data comprised information on daily fishing locations for 0.1 • spatial grids, including the fishing date, longitude and latitude, working hours (or soak time), total catch (catch weight in kg), and catch species. The life history strategies of fish differ by species, as do their environmental and habitat preferences. We used the Fish Database of Taiwan (http://fishdb.sinica.edu.tw/ (accessed on 5 January 2021)) and the FishBase database (http://www.fishbase.org (accessed on 5 January 2021)) to classify catch species as pelagic (living between the surface and bottom of the ocean) or benthic (comprising reef, demersal, and benthopelagic species).

Environmental Data and Climate Index
Daily SST data were extracted from images from the Advanced Very-High-Resolution Radiometer (AVHRR) cross-track scanning system developed by the US National Oceanic and Atmospheric Administration (NOAA). The AVHRR has a spatial resolution of 1.1 km, and the high-resolution picture transmissions are received at a ground station at National Ocean Taiwan University. Net primary production (NPP) data were estimated using the vertically generalized production model (VGPM) from the ocean productivity website of Oregon State University; the website is based on data from the Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua, which has a spatial resolution of 9 km.
The ENSO, Western Pacific Oscillation (WPO), and Pacific Decadal Oscillation (PDO) are large-scale interannual and decadal climate variabilities that develop from air-sea interactions in the Pacific Ocean [38]. In this study, the ENSO was represented by the Oceanic Niño Index, which is based on a 3-month running mean of SST anomalies in the Niño 3.4 region (5 • N-5 • S, 120 • W-170 • W) between 1971 and 2000. We used the WPO index to represent the WPO event. Based on data from the NOAA Earth System Research Laboratories (ESRL; https://www.esrl.noaa.gov/psd/data/climateindices/list/ (accessed on 5 January 2021)), the WPO index is an indicator of changes in sea surface pressure in the western Pacific Ocean in winter; such changes affect the jet stream in the western and Central Pacific, altering water temperature and rainfall patterns [39]. The PDO was represented by the PDO index, the first empirical orthogonal function characterizing detrended SST anomalies in the Pacific north of 20 • N [39]; PDO data were downloaded from the NOAA ESRL website.

Upwelling Size
The size of and variations in upwelling were estimated by examining the differences in water temperature between the upwelling and the non-upwelling regions. Upwelling regions were defined as regions with SSTs 1 • C to 3 • C lower than the monthly mean SST in the surrounding waters [1]. In images, each pixel represented 1.1 km 2 , and upwelling size was determined in accordance with the following equations used by Santos et al. [40]: where U I SST(lon,lat) is the upwelling strength of each longitude and latitude measured, SST (lon,lat) is the SST at each longitude and latitude ( • C), and SST mean is the mean SST of the study area ( • C).

PPR and Contribution to PPR
To understand the PPR, characteristics of marine environments, and fishery resource structure in the waters around the TB, we gridded the monthly environmental and fishery data onto eighty 0.2 • spatial grids. All 80 grids were classified by habitat type (upwelling, frontal, or non-upwelling) according to results reported by Lan et al. [6], who detected summertime upwelling and thermal front locations using long-term SST data ( Figure 1).
PPR was used to reflect fishery catch in comparable energy units (t C·month −1 ) for the 0.2 • spatial grids. PPR is the product of the carbon-converted catch weight and the conversion ratio for the trophic level of the taxa involved [36,41]: where C i is the catch weight (ton), TE is the trophic transfer efficiency, and TL i is the trophic level of species i. The conversion of wet weight to carbon for the calculation of trophic transfer efficiency was based on a conservative 9:1 ratio [33,42]. Data of trophic level (TL) were extracted from FishBase.
Contribution to PPR is the percentage in PPR attributable to a single fish species. This contribution for a given year is calculated as follows: where cPPR i is the percentage contribution of catch species i to overall PPR in year j, and PPR ij is the PPR of species i in year j.

Statistical Models of Fishery Resources in Each Habitat
Generalized additive models (GAMs) were used to express the effects of environmental and climatic factors on the variations in the PPR of pelagic and benthic species in the three TB habitat types. The GAMs were constructed using the GAM function of the mgcv package in R [43]. A semiparametric approach to analyzing multiple nonlinear relationships between covariates and response variables, GAMs can also explain the variance in response variables [44]. In the present study, PPR was a response variable, and environmental (SST, NPP, and upwelling size) and climatic factors (multivariate ENSO index (MEI), WPO, and PDO) were the main independent variables (p < 0.05). The GAMs were constructed as follows: where the constant c, which accounts for 10% of the mean of the overall observed PPR, was added to all PPR values because the log function cannot accept zeros; a0 is the model constant; s(x n ) is a spline smoothing function for each model covariate x n . All covariates were continuous, and the number of effective degrees of freedom was estimated for each main factor. The optimal fitting model was selected through a stepwise procedure based on minimizing the Akaike information criterion (AIC). The residual distribution was used to evaluate model fit.

Environment and Structure of Fishery Resources in TB Region
The spatial distribution of the monthly average SST and NPP in the TB region reveals that the areas of low SST and high NPP in summer closely overlap the upwelling grids, forming a banana shape ( Figure 3). Overall, the SST was lower, and the upwelling size and NPP were higher in June and July than in August and September (Figure 4a). The greatest monthly variations in environmental factors occurred in the upwelling habitat, where the minimal SST ranged from 5.64 • C to 27.61 • C, maximal NPP ranged from 926.03 to 1273.86 mg·C·m −2 ·day −1 , and upwelling area ranged from 136 to 194 km 2 per grid ( Figure 4b). The lowest SST and highest NPP were observed in frontal habitats in June also, but the upwelling spanned only 76 to 94 km 2 per grid (Figure 4c). The non-upwelling habitats (Figure 4d) had the lowest NPP and upwelling size (6-20 km 2 per grid).
Data on 311 TB species, which are caught through various fishing methods, were recorded for the study period. Pelagic and benthic species, respectively, accounted for 67.09% and 28.99% of the total catch and 77.89% and 20.46% of the PPR (Table 1). Pelagic Scombridae species constituted the highest percentage of the total catch (30.90%) and PPR (37.00%). High cPPRs were also observed for Scomber australasicus (8.03%) and Auxis rochei (7.54%), both pelagic species (Table 1). The benthic species Mene maculata accounted for 9.90% of the total catch, but the cPPRs were 3.34%. None of each benthic species was higher than 5%.

Variations in Distribution of Fishery Resources and PPR
The PPR distributions of the pelagic and benthic species in the TB region differed during the study period ( Figure 5). The PPR of pelagic species was the highest in the upwelling habitat and the northern part of the frontal habitat, as well as in the southeastern part of the non-upwelling habitat (Figure 5a). Benthic species with high PPR were observed mainly in the upwelling and frontal habitats and were both less abundant and more sparsely distributed in the non-upwelling habitat to the south (Figure 5b). Several species from the three habitats were selected for further comparison of the fishery resource structures and cPPR in the TB region. The pelagic representatives selected, Etrumeus micropus, Sardinella spp., Trachurus japonicus, Katsuwonus pelamis, Thunnus albacares, and Scombridae family. The selected benthic species are members of the Loliginidae family, as well as M. maculata and Trichiurus spp. The pelagic species from the Scombridae family and E. micropus had the highest cPPRs in the upwelling habitat, with successively lower contributions to PPR in the frontal and non-upwelling habitats ( Figure 6). The cPPRs of T. japonicus, K. pelamis, and T. albacares were higher in the frontal and non-upwelling habitats. Benthic species from the Loliginidae family had higher cPPR in the upwelling and frontal habitats than in the non-upwelling habitat, and the cPPR of M. maculata was high (5.81%) in the non-upwelling habitat.

Effects of Environmental and Climatic Factors on PPR
The data used to select the benthic and pelagic species for GAM analysis are presented in Tables 2 and 3. The maximal deviance explained of the selected benthic and pelagic species PPR was 28.5% and 22% in the frontal habitat (Table 2), 19.2% and 16.9% in the upwelling habitat, and 8.9% and 11.5% in the non-upwelling habitat. The NPP accounted for considerable portions of the variation in benthic species PPR in all three habitats, as well as that in the pelagic species in the frontal habitat (Table 3). Both SST and upwelling size were critical factors, explaining 6.41% and 6.83% of the deviance of the pelagic species PPR in the frontal habitat. Overall, the climatic factors of the MEI and PDO index were most strongly associated with the pelagic species PPR in the upwelling habitat (MEI: 10.10%) and the non-upwelling habitat (MEI: 4.74%, PDO: 5.26%).  Table 3. The deviance explained by each factor in each habitat GAM model.

s(SST) s(NPP) s(Upwelling_Size) s(WPO) s(MEI) s(PDO)
Upwelling A positive association was observed between the NPP and PPR of benthic and pelagic species in all habitats for NPP concentrations of 500-1500 mg·C·m −2 ·day −1 in all three habitats (Figures 7 and 8). Monthly NPP peaked in June and decreased gradually through September for both species types in all three habitats, whereas monthly PPR increased gradually from June to August. These results indicate that the relationship between PPR and NPP may be time-lagged ( Figure 9). Moreover, upwelling size and SST (27-29 • C) were positively associated with the benthic PPR in the frontal habitat (Figure 7c,d).
As shown in Figure 8, the MEI and PDO index were negatively associated with pelagic PPR in the upwelling and non-upwelling habitats, which revealed the differences in the characteristics of changes between years. The distribution of and differences in the annual summertime averages of the SST, NPP, and PPR of pelagic and benthic species during the negative (La Niña) and positive (El Niño) MEI events are presented in Figure 10. The annual summertime average time series indicate that pelagic PPR was highest in 2011 in La Niña year with higher NPP and relatively low in 2016 in El Niño year (Figure 10e). Otherwise, the relationships between benthic PPR and ENSO events were not clearly on time series patterns (Figure 10f).   On the other hand, the difference of annual spatial distribution reveals that higher average PPRs were noted of pelagic species in the northwest and southeast parts of the frontal habitat, northern parts of the upwelling habitat, and southern parts of the nonupwelling habitat (Figure 10c), with lower SST (Figure 10a) and higher NPP (Figure 10b) in La Niña years. Variations in the spatial distribution of the PPR of benthic species were not associated with environmental factors during ENSO events (Figure 10d).

Environmental Characteristics of Three TB Habitats
The monthly average data indicated that upwelling size and NPP were higher in June and July, and that in the waters around the TB, the SST, which is lower in summer, generally differs among the three habitats. For the summertime average SST difference between the upwelling habitat and fronts, non-upwelling is close to 0.7 and 1.5 • C, respectively. The SST was around 26 • C in TB upwelling, which was also within the SST ranges (24-26 • C) of upwelling areas shown in previous studies [1,3,6]. Furthermore, the SST gradient magnitudes of thermal front habitat calculated by Lan et al [6] were 0.1 to 0.2 • C/km and gradually disappeared in September. With other upwelling near TS, the average SST differences between upwelling and non-upwelling areas were 2 to 4 • C in the East China Sea and Dong Sheng upwelling [1,45].
Compared to Tang et al. [1], our results showing the difference of SST between upwelling and non-upwelling areas in TB presented the same patterns, which were approximately 1-3 • C, averaging 2.35 • C [1,46]. In September, as the SCSC and southwest monsoon weaken, the northeast monsoon and CCC begin to strengthen [6,47]. Lin et al. [48] discovered that rather than being closely correlated with variations in upwelling strength, wind direction and speed were mainly associated with topography, tides, and ocean currents [3,46,49]. The southwest monsoon drives the SCSC and KBC into the TS, causing the velocity of surface currents to exceed that of bottom currents [47,50]. SST increases gradually from upwelling areas to non-upwelling areas, and frontal habitats are affected by the upwelling of cold bottom water and downwelling of warm surface water [6,46,51], an effect that was strongest in July and weakest in September in this study.
The NPP, which also had different concentrations in three habitats, and NPP difference compared with upwelling habitats were closed to 238 and 590 mg·C·m −2 ·day −1 in thermal front and non-upwelling habitats, respectively. The monthly mean NPPs of MODIS Aqua calculated by VGPM were in the ranges of 500-1500 mg C m -2 day -1 in TS [52]. The NPP ranges of China coastal upwelling is 600-2500 mg C m -2 day -1 ; in Taiwan, northeastern upwelling is 500-1000 mg C m -2 day -1 [52,53]. Rather SST or NPP in TB and other upwellings around TS present the obvious difference between upwelling and non-upwelling habitat. We suggest which could generalize into SST deviation at a lowest of 1 • C and NPP deviation at a lowest of 500 mg C m -2 day -1 . These differences not only indicate the geography status but also show the impact of upwelling mechanism. The northeastern upwelling in TS was mainly caused by shelf break and wind stress curl; China coastal upwelling along TS was caused by alongshore winds, and TB upwelling was suggested to be caused by shelf break and background flow [3,7,23,45].
The NPP data indicated low nutrient concentrations in the non-upwelling area, whereas abundant nutrients such as phosphates and nitrates were present in the water upwelled from lower to upper layers or even to the surface layer [45,47]. Spatial distribution data suggest that the high PP, attributable to upwelling events and the accompanying nutritional richness, was responsible for the negative correlation between SST and NPP. In addition, in the first comparison of the satellite-derived VGPM estimates with in situ PP in the TS, Lan et al. [52] reported high correlation and low root-mean-square deviation. Satellite SST represents temperature only in the uppermost ocean layer and P B opt , being a function of SST, and differs in the lower layers, providing different photosynthetic efficiency [27]. The VGPM is one of the most widely known and applied depth-integrated/wavelength-integrated models. However, it has rarely been applied to coastal waters [54]. Thus, improvements in the accuracy of Chl-a from other ocean color sensors, including Medium Spectral Resolution Imaging Spectrometer and Ocean Color Climate Change Initiative data, will ultimately lead to an improvement in satellite PP algorithms for further research.

Fishery Resource Structure and PPR in Three TB Habitats
Both the fishing methods and fishery resource structure in the TB region are complex and variegated. The TS is located on the northern boundary of the Kuroshio triangle; it is characterized by high diversity of corals and reef species, and the currents and water depth create a suitable habitat for migratory species [55][56][57][58]. The distribution of top predators is dependent on that of anchovies and other forage species, which cluster in habitats located in thermal fronts [17]. Tiedemann and Brehmer [59] observed that larval and juvenile fish are divided into groups living in upwellings or thermal fronts. We noted that both fishery resource structure and PPR distributions in the TB differed by habitat.
In the present study, the upwelling habitat had high NPP and was dominated by cold-water pelagic species (i.e., E. micropus and Scombridae species) and benthic species (i.e., Loliginidae species). Attracted by the abundance of forage species, Scombridae species inhabit low-temperature areas with high NPP [60,61]. The main higher trophic level species (TL > 4.0) were Trichiurus spp. (TL = 4.42), followed by Loliginidae (TL = 4.0), Scombridae (TL = 4.0) and Etrumeus micropus (TL = 3.5). The Trichiurus spp. and Scombridae in TS main prey the Cephalopoda, which also prey small reef/demersal fish and small pelagic fish, respectively [62][63][64]. The species variations between three habitats were not only effected by the environment factors but also connected to the diet and feeding habits.
Furthermore, the PPR contributions of high trophic level species of Trichiurus spp. and Katsuwonus pelamis (TL = 4.03) were increased in the thermal front habitats, followed by Loliginidae, Etrumeus micropus and Trachurus japonicas (TL = 3.4). Tanabe [65] shows that the Trichiurus spp. and Loliginidae were the two main forage species of Katsuwonus pelamis. Trachurus japonicas is a small pelagic fish, and its diet is partially similar to Scombridae [66]. In the non-upwelling habitat, Thunnus albacares and Mene maculate express the higher composition, suggesting co-variation with the predator-prey relationships [67].
In addition, some tropical species such as T. albacares, Coryphaena hippurus, and K. pelamis appeared in the non-upwelling habitat. Chen [60] noted that the summer southwest monsoon causes a gradual rise in SST, and that upwelling periods characterized by high NPP correspond to the periods of the highest T. albacares and K. pelamis catches in non-upwelling habitats. In the present study, the fishery resource structure in the frontal habitat indicated that it is a transition zone between upwelling and non-upwelling habitats, a premise supported by the presence of both warm-and cold-water species, including E. micropus, K. pelamis, Scomberomorus commerson, T. japonicus, and Loliginidae species. Sabatini and Martos [17] and Sassa et al. [68] have reported that T. japonicus and Loliginidae species concentrate along the edge of thermal fronts. Wang et al. [69] found that frontal areas along northeastern Taiwan are dominated by immature fish from the Loliginidae family, whereas their adult counterparts appear in coastal waters and upwelling areas.

Explanatory Variables in GAM Analysis
The GAM analysis revealed that NPP accounted for a considerable portion of the variation in the benthic species in all three habitats, as well as that in the pelagic species in the frontal habitat. Moreover, the variation in PPR with NPP in the upwelling and frontal habitats may be time-lagged by 1 to 2 months. Higher NPP increases zooplankton blooms, directly increasing the first-level consumer resources [70,71]. Furthermore, the downstream secondary production creates an attractive habitat for species, such as those of the Scombridae and Loliginidae families and T. albacares, which occupy high trophic levels [28,53]. Therefore, the effect of high chlorophyll-a concentrations on NPP appears to be time-lagged because of trophic transfer efficiency [32,52,72,73]. Furthermore, the cPPR of benthic species was affected by variations in NPP. Studies have suggested that high spatiotemporal variability in NPP is positively correlated with demersal production [74,75]. The benthic PPR was not correlated with climate index but was highly positively correlated with NPP, indicating that its effects on NPP and upwelling events are asynchronous. Notably, the effects of dominant ENSO events on fishery catch were inconsistent. Furthermore, no relationships between the ENSO and NPP were noted. Tzeng et al. [76] observed that not all El Niño/La Niña events are similarly correlated with NPP in the eastern TS, stating that other factors may influence the spatial distribution of NPP. Wu et al. [77] inferred that the ENSO does not exert one single effect on summertime upwelling in the TS, and that the interaction between the East Asian monsoon systems in the China seas may play a role [78,79].
Our results indicate that pelagic PPR in the TB region is primarily affected by environmental and climatic factors. Song et al. [80] indicated that the SST is the factor most contributing to spawning-feeding migration in the TS. Ma et al. [81] reported that small pelagic species, including the chub mackerel and Pacific sardine, favor the low temperatures of the East China Sea. Climate variabilities have changed the marine environment in the TS; specifically, ENSO events influence the intensity of the upwelling in the TS [82]. The weaker wind field during El Niño events explains why summertime SST is higher during those events than during La Niña [12,25]. In essence, the changes in SST and upwelling intensity caused by climatic factors in the TB are the main factors affecting the migration patterns of the pelagic species in the three habitats.

Effects of Climate Change on the TB
The thermal fronts in the TB are characterized by two frontal belts on the northern and southern edges. Studies have indicated that the formation of the thermal front along the southern edge is linked to the interaction between topography and upwelling phenomena, whereas the thermal front along the northern edge may be tidally generated [51,83]. We compared the variations in the cPPR of pelagic and benthic species during El Niño and La Niña events (Table 4). During El Niño events, catches of Scombridae species decreased (Table 4), whereas catches of low-trophic-level species, such as E. micropus and Sardinella spp., increased. Moreover, pelagic PPR exhibited a trend of decline during these events ( Figure 10). The cPPR of the pelagic species was higher during La Niña events, especially in the southern and northern part of the non-upwelling and upwelling habitats, respectively ( Figure 10c). Moreover, cold-water species (Scombridae and Carangidae species and E. micropus) increased in abundance, as did warm-water species. Annual average time series obviously showed that the pelagic PPR was higher at La Nina (Figure 10e). The stronger El Niño event from 1997 to 1998 changed the dynamics of the summer monsoon and the mechanism of upwelling in the TS. Furthermore, it affected the composition and abundance of high-trophic-level pelagic species [25,78]. The PDO has been shown to affect ocean productivity and eddy currents in the highand mid-latitude regions of the North Pacific [76,84]. These effects are in turn associated with changes in marine ecosystems and fisheries [20,85,86]. The PDO index was a crucial factor in our GAM results, particularly those for the pelagic species in the upwelling and non-upwelling habitats. The WPO has been defined as a regional climate mode or subseasonal mode in which teleconnection patterns are expressed in a common region that has e-folding timescales between 7 and 10 days [87,88]. However, how the WPO alters climate patterns in the TS remains relatively unexplored. Baxter and Nigam [89] speculated that the WPO exerts stronger effects on air-sea flux in air pressure and rainfall. Thus, collecting more long-term fishery data and expanding the time scale on which the impacts of the PDO are investigated (i.e., beyond the decadal scale) are warranted in future research.

Conclusions
In the present analysis of fishery resource structure and the spatial distribution of PPR in upwelling, frontal, and non-upwelling habitats in the TB region in summer, variations in benthic PPR were correlated with changes in NPP. Members of the Loliginidae family, as well as M. maculata and Trichiurus spp., were the major benthic species in all three habitats. Variations in the pelagic PPR were influenced by ENSO events, and increased upwelling size contributed to lower SST and higher NPP in the TB region during La Niña events. High cPPR during La Niña events was observed for the pelagic species in the upwelling habitat, mainly those in the families of Scombridae (including A. rochei) and Carangidae. By contrast, during El Niño events, the dominant species were E. micropus and S. commerson. The frontal habitat marked the zone of transition between the two other habitats, and members of the Scombridae and Carangidae families and E. micropus had high cPPR in La Niña years. The high cPPR species during El Niño years were T. japonicus and K. pelamis. During La Niña events, the dominant pelagic species in the non-upwelling habitat were the same as those in the frontal habitat. During El Niño events, A. rochei rochei was the dominant species ( Figure 11).
In the future, we hope to evaluate the performance of our models through various habitat classification methods. The global trend of sea temperature rise exerts substantial impacts on fisheries and ecosystems [77,90]. Shifts in environmental characteristics under the impacts of climate change mean that some habitats are no longer suitable for certain species; such impacts affect the entire ecosystem. The present research period spanned