The Visayan Sea is located in the central Philippine archipelago (11.4328° N, 123.5280° E) with a total surface area of approximately 12,521.8 km² and 16,728.5 km² if municipal waters outside the Visayan Sea domain are included, respectively ( Figure 1 ). For the purpose of this study, the Visayan Sea baseline Ecopath model includes the delineated commercial water area and municipal waters (both domain and non-domain) of all 33 coastal municipalities and cities that surround the Visayan Sea. Due to a lack of information on the spatial distributional patterns of fishing effort, fisheries activities are hereby initially assumed to occur throughout the entire extent of the delineated model area.

Due to its central location and physical oceanographic attributes, the Visayan Sea has long been one of the country’s major fishing grounds that caters to both the commercial and municipal fisheries sectors (as classified under the Philippine Fisheries Code of 1998). Historically, the commercial fishing sector in the Visayan Sea primarily includes the purse seine, ring net, bag nets, Danish seines, and trawling operations that involve the use of vessels with a capacity > 3 gross tonnage (G.T.). On the other hand, the municipal sector includes the artisanal fisheries and fishing operations using vessels with a capacity< 3 G.T. In the mid-1990s, the Visayan Sea was already the highest contributor to the country’s overall commercial fisheries landings (at 120 thousand metric tons/year) and ranked the 3rd highest contributor to the country’s municipal fisheries sector (at 89 thousand metric tons/year) (DA-BFAR, 2000).

Surveys conducted by the National Stock Assessment Program (NSAP) between 1998 and 2002 identified a total of 96 families of aquatic fauna (90 vertebrates and 6 invertebrates) in the fisheries landings from the Visayan Sea (Guanco et al., 2009). These can be classified further into 149 genera (138 vertebrates and 11 invertebrates) and 282 species (264 vertebrates and 18 invertebrates). Fisheries monitoring data from the NSAP in the past decade show that the sardines (Clupeidae) and mackerels (Scombridae) comprise 50% of the total small pelagic fish landings, while goatfishes (Mullidae) and monocle breams (Nemipteridae) comprise about 70% of the soft-bottom demersal fish landings. At the same time, the landings of invertebrates such as tropical squids (Loliginidae) and of blue swimming crabs (Portunidae) have likewise grown steadily from 4% to 12% of the total landings.

For this model, a total of 33 groups were identified, including 30 consumers, two primary producers and one detritus group ( Table 1 ). The groupings were initially based on the documented fish and invertebrate species in the Visayan Sea (Armada, 2004; Guanco et al., 2009; Del Norte-Campos et al., 2019) and further defined based on major taxonomic grouping (i.e., fish, invertebrate), ecological function (i.e., primary producer, consumer), and similarities in food types and feeding behavior (i.e., predator/hunter, grazer/herbivore, omnivore, piscivore, planktivore) classifications in Fishbase (Froese and Pauly, 2020). In addition, they were further classified based on a general knowledge of their habitat associations (i.e., pelagic, soft-bottom/demersal, reef-associated) in the Visayan Sea as informed by the fisheries landings surveys (Mesa, 2019), fisheries-independent trawl surveys (Ampoyos-Arinque, 2018), and coral reef monitoring surveys (Campos et al., 2020c; Campos et al., 2020a; Campos et al., 2020d; Campos et al., 2020b). Further, fish and invertebrate taxa of known major fishery importance (e.g., small pelagics, demersal fish, rabbitfish, marine crabs, and squids) comprise distinct functional groups. In addition, pelagic predator fishes (primarily Scomberomorous species) were split into juvenile and adult stanzas to incorporate ontogenetic shifts in diets and other life history parameters. The three small pelagic fish groups, namely the sardines, mackerels, and scads, were likewise split into juvenile and adult stanzas in anticipation of future explorations of the potential impacts of existing fisheries policies ¹ , which may have differential impacts on the juvenile and adult stages.

The types and numbers of fishing gears and operations in the Visayan Sea are described in Armada (2004) and updated in Bacalso (2019). They include several variations of commercial and municipal fishing gears and operations such as gillnets, hook and lines, jigs, pots and traps, seines and dragnets, spears, impounding nets, falling nets, lift nets, surrounding nets, and even unconventional fishing methods such as blast fishing. For this model, they are classified into 34 gear groupings or “Ecofleets” according to basic gear type, function, target species, and/or relevance to existing national fisheries policies and local regulations (e.g., use of fine mesh nets, active gears, and destructive gears) ( Table 2 ). In terms of number of units, the municipal fishing operations clearly outnumber the commercial operations. Owing to their higher catching efficiencies, however, the commercial operations contributed 56% of the estimated total catch from the Visayan Sea in 1996-1997.

A mass-balanced trophic model of the Visayan Sea was constructed using the Ecopath with Ecosim software (www.Ecopath.org) to describe the energy pathways and ecological dynamics therein. Mass-balanced system models (Christensen and Walters, 2004) assume that flows of biomass (or energy) balance around trophically-interacting ecological compartments within a defined time period. In this case, the baseline model of the Visayan Sea represents average values of fisheries and biological information from a comprehensive year-long fisheries assessment study conducted in the years 1996-1997 (Armada, 2004).

The parametrization of the model is based on trophic groups production and consumption, which can be expressed formally as two master equations (Christensen et al., 2005; Christensen et al., 2008):

where (P/B)_i is the production/biomass ratio of i which is equal to the coefficient of total mortality (Z) under steady state conditions; B_i is the biomass of prey group i; Y_i is the fishing mortality of group i; B_j is the biomass of predator j; (Q/B)_j is the consumption/biomass ratio of predator j; DC_ji is the fraction of prey i in the diet of predator j; E_i is the net migration rate of group i; BA_i is the biomass accumulation of group i; EE_i is the ecotrophic efficiency or proportion of production that goes to predation, catches or exports out of the system; and (1-EE_i ) is the other mortality (not caused by either predation or fishing).

For the second master equation:

Q_i is the consumption, P_i is the production, R_i is the respiration, and U_i is the unassimilated food of group i. For each functional group in the model, biomass (t/km²), production/biomass ratio (/year), consumption/biomass ratio (/year), diet composition (%) and contribution to the fisheries catches (t/km²/yr¹) are required inputs. Information used to construct the Visayan Sea baseline Ecopath model were sourced primarily from the 1996-1997 Fisheries Assessment Survey report (Armada, 2004) and succeeding local fisheries surveys and assessments conducted by the NSAP (Guanco et al., 2009; Mesa, 2019). Data “pedigree” scores were assigned via the built-in model Pedigree function to classify the confidence of data inputs based on their sources (i.e., site based/high precision sampling, limited/incomplete sampling, indirect/approximation methods, general knowledge, and from other models of similar systems) (see Supplement Table S1 ).

Results from the BFAR-led exploratory trawl surveys from 2003 and from earlier exploratory trawl and marine fisheries surveys in the Visayan Sea (Warfel and Manacop, 1950; Rasalan, 1957; Aprieto and Villoso, 1979) were used as basic reference for initial inputs of demersal fish biomass. Initial inputs of the pelagic fish group biomasses were either estimated indirectly using the relationship B = Y/F where F is the fishing mortality rate (/year) and Y refers to the estimated fisheries landings (t/km²/year) or is left to be estimated by Ecopath (for example, the biomass of the juvenile stanzas and most invertebrate groups). On the other hand, biomass inputs of the reef-associated species were sourced from the reef-based fish visual census surveys in the Visayan Sea, which report a range of reef fish biomass estimates from 1.2 – 26.7 t/km² (Campos, 2018). The reef habitats cover less than 10% of the delineated area for the model. However, the proportion of habitat area occupied by reef and reef-associated functional groups was assigned a relative value of 0.2 to account for movements between reef and non-reef habitats. The other benthic, soft-bottom demersal and pelagic functional groups were assumed to occur over the entire model area. Phytoplankton biomass was estimated indirectly from satellite-based chlorophyll a data (source: NOAA/NCEP Global Forecast System (GFS) Atmospheric Models) and using conversion factors in Opitz (1996). Biomass of the detritus compartment was likewise estimated empirically sensu Pauly et al. (1993b).

Production/biomass ratios for groups exploited by fisheries were weighted averages relative to harvest weight of the instantaneous rate of total mortality (Z) estimates from published and unpublished stock assessments of fish and invertebrate stocks in the Visayan Sea (Guanco, 1991; Armada, 2004; Guanco et al., 2009; Romero, 2009; Tajolosa, 2011), as well as, biological parameter estimates of the same representative fish species from other Philippine fishing grounds around the same time period (Cinco and Silvestre, 1994; Armada et al., 2003). Total mortality of fish groups not covered in fisheries assessments were initially assumed equal to the instantaneous rate of natural mortality (M) (Pauly, 1980; Djabali et al., 1993). The consumption/biomass ratios for fish groups were estimated using an empirical formula on predicting food consumption rates (Pauly et al., 1990; Palomares and Pauly, 1998). Initial inputs of Q/B for the invertebrate groups were either sourced from the literature (Pauly et al., 1993a) or left to be estimated by the model. The initial Biomass, P/B, and Q/B inputs were then adjusted accordingly during model balancing and parametrization. Diet inputs for the demersal fish groups were sourced from a quantitative diet composition study in the Visayan Sea (Mequila and Campos, 2007). For the rest of the consumer groups, published diets of the same species from similar systems (i.e., shallow bay fishery systems) in the Philippines and elsewhere in the tropics were consulted. The complete diet composition matrix is provided in the Supplementary Material (Table S2) .

The expanded baseline fisheries inventory data and gear CPUEs (kg per unit of operation) were used to estimate the total fishing effort and the total catch (t/km²) of the exploited groups by the different Ecofleets ( Supplement Table S3 ). The total estimated catches were then apportioned to the different functional groups based on reported species and percentage composition in the catch (Armada, 2004; Guanco et al., 2009).

Ecosim was used to simulate changes in biomass and catch over time with fishing effort of the different Ecofleets and several environmental parameters as drivers. Based on the initial parameters of the Ecopath master equation (Equation 1), a series of coupled differential equations are used to express the biomass dynamics in Ecosim:

where dB_i/dt is the biomass variation of group i through time; g_i is the net growth efficiency (equivalent to the P/Q ratio); I_i is the immigration rate; M_i is the natural mortality rate; F_i is the fishing mortality rate; and e_i is the emigration rate. The calculation of the consumption rates Q and the vulnerable and non-vulnerable components of the prey biomass B are based upon the “foraging arena” concept, which describes the foraging time of an animal as a tradeoff between feeding and predation risk (Ahrens et al., 2012).

In Ecosim, trophic control mechanisms between the model’s functional groups are defined by the vulnerability term v (Christensen et al., 2005; Christensen et al., 2008), which can be set between 1 and ∞. A low v setting (i.e., close to 1) indicates a bottom-up control whereby an increase in the predator biomass does not result in a substantial increase on the prey’s predation mortality. In contrast, a high v setting (e.g., 100) indicates a top-down control whereby an increase in predator biomass results in a substantial increase in the predation mortality of the prey. A high v likewise indicates that the predator biomass is far below its carrying capacity (Guénette et al., 2008). The default v is 2, which indicates neither a top-down nor bottom-up control (or a mixed control setting).

The basic input and output parameters of the balanced Visayan Sea baseline Ecopath model are presented in Table 3 . Model PRE- and POST-BAL diagnostics (see Supplement Figures S1–S4 ) for the Biomass, P/B, and Q/B input and output parameters improved (regression fitting, R²), while the P/Q estimates were likewise constrained under 0.3 after model balancing and parametrization. Overall, the final Visayan Sea baseline model adhere to the basic rules and principles in system ecology. The EE values for most invertebrate groups, the small pelagics (particularly the juvenile stanzas), and the small demersal fishes are high (close to 1) indicating that these groups are highly utilized in the system as prey for the different consumers, as targets of the fishery, or both. The model estimated F/Z (exploitation ratios) for the small pelagic groups (0.5-0.67) and small demersal fishes (0.51) are also within the ranges reported in single-species assessments from the Visayan Sea and similarly highly exploited Philippine fishing grounds (Ingles and Pauly, 1984).

The flow diagram ( Figure 3 ) visualizes the trophic relationships and relative biomass proportions of the different groups in the model. The energy from TL 1 is first utilized by several planktivores, herbivores, and detritivores before reaching the higher order consumers. In terms of energy utilization, there seems to be no clear delineation between the pelagic and demersal groups since demersal prey (such as shrimps and other benthic invertebrates) are also preyed upon by pelagic consumers. At the same time, zooplankton and pelagic prey fish groups (e.g., anchovies and small pelagic fish, especially the juvenile stanzas) are also consumed by squids and the demersal predators. The shallow nature of the Visayan Sea is seen as the main physical characteristic that facilitates this interaction.

Selected outputs from the system summary statistics of the Visayan Sea baseline Ecopath model are shown in Table 4 . System size in terms of total system throughput (TST) is in the mid-range compared to the throughput of the 3 other published trophic models of Philippine fishing grounds that represent system states within 5 years from the Visayan Sea baseline model. The Visayan Sea has the highest total catch (t/km²/yr) followed closely by San Miguel Bay (Bundy and Pauly, 2001). The mean trophic level of the fishery (TL_C) is computed as the weighted average TL of harvested groups, and it likewise serves as an index of fishing impact (Christensen et al., 2005; Pauly et al., 1998). The TL_C for the Visayan Sea is very similar to the reported TL_C of San Miguel Bay and Sapian Bay (Bacalso and Armada, 2004), which could indicate similarities in the overall composition of the catches. Meanwhile, the gross efficiency of the catch is computed as the sum of all realized fisheries catches relative to the total primary production (catch/net p.p.), and it tends to be higher in systems targeting low trophic level resources (Christensen et al., 2005). The value for the Visayan Sea comes second to San Miguel Bay, followed by Sapian Bay, and finally, San Pedro Bay. Based on these two indices alone – TL_C and catch/net p.p. – we can say that the Visayan Sea, San Miguel Bay and Sapian Bay were similar in terms of fisheries exploitation during the 1990s. This finding likewise validates the contention that San Pedro Bay was not as heavily exploited as the other traditional fishing grounds in the country during that time (Campos, 2003).

The total primary production/total biomass ratio (PP/B) is an indicator of system maturity (Christensen et al., 2005). The ratio is expected to decline as biomass accumulates over time or as a system matures. The PP/B ratio for the Visayan Sea is closest to San Miguel Bay and Sapian Bay, which could again indicate their similarities in terms of the types and rates of fisheries extraction. The total biomass/total throughput ratio (B/T) is another indicator of system maturity (Christensen, 1995) and increases as a system matures. The values of all 4 fishing grounds hereby compared are similar with each other. Lastly, the Finn’s Cycling Index (FCI) represents the fraction of an ecosystem’s throughput that is recycled. It is an indicator of a system’s ability to maintain structure and integrity (Ulanowicz, 1984). The Visayan Sea has the highest FCI compared to the other models, which indicates the Visayan Sea’s higher resilience toward disturbance compared to the other systems.

The computed sum of squares and AICc of all the models that were generated under different hypotheses are shown in Figure 4 . The application of fishing effort time-series alone (H5) results in the largest disagreement between the historical catch and biomass with the model simulations. As likewise shown in Figures 5 , 6 , even the best-fitting model under this scenario usually either underestimates or overestimates the catch and biomass of several groups. It further predicted a collapse of the biomass of mackerels and scads in the last years of the simulation ( Figures 5B, C ) when the fishing effort of the commercial ring nets and most municipal fishing gears were at their peak ( Figure 2 ). It also predicted a corresponding crash in these groups’ catches ( Figures 6B, C ) as well as in the catch of squids ( Figure 6E ).

On the other hand, the models with both fishing effort time-series drivers and vulnerability setting adjustments (H7 to H12) outperformed the rest by producing lower values of model SS and AICc overall ( Figure 4 ). The best-fitting model was produced when time-series of sea surface temperature (sst) was added as a forcing function to the consumer search rates along with a search for primary productivity anomalies (H11) ( Figure 4 ). This resulted in the largest decrease of model SS and AICc from the baseline scenario ( Table 5 ). This model also achieved the best fit for the biomass of mackerels ( Figure 5B ) and the biomass and catch of crabs ( Figures 5F , 6F ). The best fit for the catches of sardines and mackerels were achieved when time-series of rainfall was applied along with a search for primary productivity anomalies ( Figures 6A, B ). However, no model was able to reduce the SS for the sardine biomass below the baseline scenario. In other words, neither fishing nor any of the environmental drivers was able to better explain statistically the historical biomass trend of the sardines, which appear to be stable over time ( Figure 5A ). All alternative models significantly improved the fits for the small demersal fish biomass. However, the best fit was achieved when time-series of chlorophyll-a was applied as a forcing function along with a search for primary productivity anomalies ( Figure 5D ). Similarly, best fits for the catch of this group ( Figure 6D ) were obtained when chlorophyll-a was applied as forcing function. For the scads, no model was able to reduce the biomass and catch sum of squares below the baseline scenario. Nonetheless, better fits were obtained visually with the observed catches when fishing effort drivers, vulnerability setting adjustments and PP anomaly search were applied ( Figure 6C ). All the alternative models with fishing effort drivers, vulnerability adjustments, PP anomaly search, and environmental forcing functions improved the fits to the historical squids biomass. However, fishing effort as the sole driver resulted in a model with the largest % reduction in SS ( Table 5 ). On the other hand, no alternative model was able to reduce the catch SS for this group compared to the baseline scenario. Visually, however, better fits were obtained when time-series of SST was applied as a forcing function ( Figure 6E ).

Overall, the results show that majority of the groups are highly sensitive to the effects of variable fishing effort, but trophic interactions play an important role in modulating these effects. The adjusted vulnerability settings from the best-fitting model show the most principal trophic control mechanisms in the Visayan Sea ( Table 6 ). The adult stanza of the sardines, the juvenile stanza of the mackerels, the adult stanza of the scads, and the squids are shown to be bottom-up controlled by their prey. On the other hand, the adult mackerels, the juvenile stanza of the scads, and the reef-associated predator group 1 (jacks and barracudas) exert a top-down control on their prey.

Local environmental variabilities also influence the biomass and catch of the different groups, with several groups being more sensitive to a particular parameter than the others. Interestingly, the influence of primary productivity anomalies is not exclusive to the plankton-feeding small pelagic groups as shown by its contributions to improving the individual fits to biomass and catch time-series of several groups. The primary productivity anomalies from the best-fitting model overall were shown to be positively correlated with the Pacific Decadal Oscillation (PDO) and negatively correlated with the Southern Oscillation Index (SOI) ( Figure 7 ), indicating a possible bottom-up influence of large-scale and longer-term climatic variabilities in the Visayan Sea ecosystem. Among the local environmental parameters, only rainfall showed a weak but nonetheless significant negative correlation with the primary productivity anomalies ( Figure 7 ). Subsequently, rainfall is negatively correlated and positively correlated with PDO and SOI, respectively.

Applying the ecosystem modeling approach allowed us to assess the status of the Visayan Sea as a whole in relation to the influence of both fisheries and environmental parameters on the system’s trophic interactions. From the model outputs of ecosystem network indices, we were able to identify system-level expressions of biomass extraction by the fishery (e.g., trophic level of the catch and G.E. of the fishery) and make inferences on the Visayan Sea’s status and development (Odum, 1969; Christensen, 1995). Importantly, these indices allowed us to compare the status of the Visayan Sea with other Philippine fishing grounds, and showed that despite their differences in geographic location and size, systems with similar fisheries extraction patterns produce closely similar indices. As shown by Heymans et al. (2014), marine ecosystems of the same type generally share similarities in system and network indices, but may either vary or correspond with each other with changing fishing intensity. This approach is very useful as it provides the opportunity to compare not only discrete ecosystems (Del Solar et al., 2022), but also the state of the same ecosystem across time with changes in fishing exploitation patterns (Neira et al., 2004) and under different climate regimes (Taylor et al., 2008).

Our results and those from the other Ecopath models of Philippine fishing grounds compared herewith confirm the considerably high impact of fishing activities even from two decades ago. They likewise reflect similar network indices and mortality rates reported in other highly exploited fishing grounds elsewhere in tropical Asia (e.g., Liew and Chan, 1987; Christensen, 1998; Chen et al., 2015). Going further back in time to the earliest exploratory trawl surveys in the Visayan Sea (Aprieto and Villoso, 1979), the demersal standing stock was estimated at around three times this baseline model’s estimate. In addition, what was considered as “trash fish” in the catch composition then (e.g., Apogonidae, Centriscidae, Centropomidae, and Pomacentridae) (Ordonez, 1985) feature more prominently in the catch composition of the trawl surveys from the last decade (Ampoyos-Arinque, 2018). Based on these findings, we can infer that the Visayan Sea ecosystem has undergone changes long before the first comprehensive fisheries assessment was conducted in the area (late 1990s). While we cannot determine the specific cause of these past changes, they have often been attributed almost exclusively to overfishing, as is the case with several other Philippine fishing grounds that have shown stock and catch declines over the past decades (Green et al., 2003; Barut et al., 2004; Anticamara and Go, 2016).

Lately, however, there is a growing recognition of the effects of climate and environmental variability on stock abundance in the Philippines (Damatac and Santos, 2016; Geronimo, 2018) and on the fisheries in general (Gopal and Anbumozhi, 2019; Suh and Pomeroy, 2020; Macusi et al., 2021). In a recent study, inter- and intra-annual environmental variabilities were shown to have significant influence on the species and gear CPUEs in the Visayan Sea (Bacalso et al., 2023). Through the ecosystem modeling approach, the present study was able to demonstrate the discrete and combined effects of fishing and the environment on the catches and biomass of the important fisheries groups in the Visayan Sea. The model fitting to time series under different hypotheses has shown that fishing effort trends in the Visayan Sea cause the largest deviation from the baseline model scenario ( Figure 4 ), which suggests that indeed, the fishery has considerable influence on the historical catch and biomass. However, unrealistic trends were predicted with fishing effort as the sole driver of the system, such as scenarios of biomass and catch collapse in the later years when the effort of most fishing gears were at their peak. The consideration of trophic interactions significantly improved the overall model fits, while the addition of environmental parameter time series resulted in further refinements to the fits for several of the groups’ catch and biomass ( Figures 5 , 6 ).

Consistent with the results shown in Bacalso et al. (2023), the present study showed the significant influence of SST and primary productivity as the parameters that produced the best-fitting model overall. Likewise, the present study confirmed the important and seemingly favorable influence of SST on marine crabs in the Visayan Sea, particularly the blue swimming crabs whose species CPUEs show concurrence with the SST trend (Bacalso et al., 2023) due to a combination of biological and behavioral responses to increasing temperatures (De Lestang et al., 2003a; De Lestang et al., 2003b; Johnston et al., 2021). In the present study, this increase also appears to coincide with an increase in the effort of crab fishing gears and in the overall fishing effort, particularly from the municipal fisheries sector from 2012 onwards ( Figure 2 ). However, the catch of marine crabs appears to have already reached a peak in 2015 followed by progressively decreasing catches in the succeeding years ( Figure 6F ) despite the high fishing effort. This is an indication that the CPUE increase on the marine crabs was not sustainable in the long term. In contrast to the increasing overall trend of the municipal fishing effort, the effort of most commercial fishing gears showed steady declines since the passing of the Philippine Fisheries Code ² in 1998 that set regulations and area limitations on the commercial fisheries activities. Further, a moratorium ³ on the issuance of commercial fishing vessel and gear licenses was decreed in 2014, which was followed by tighter penalties for fisheries violations in the amended Fisheries Code ⁴ in 2015. In addition, a country-wide ban on Danish seine operations ⁵ was declared in 2013, which apparently resulted in further declines in their fishing effort. Coinciding with the increasing SST trend and the generally contrasting effort trends of the two fisheries sectors are the declining trends in the catch and biomass of mackerels, small demersal fishes and squids. Despite the absence of long-term biomass time-series for the predator groups, the model likewise predicted similar declining trends for these groups within the 20-year simulation period. The considerably high (> 100) adjusted vulnerability setting for one of the predator groups in the system – the jacks and barracudas – indicate that their initial biomass was already far below the carrying capacity, and that subsequent changes are expected to result in corresponding changes in the mortality rate of their prey (Guénette et al., 2008). In this case, the declining predator biomass resulted in a corresponding decrease in predation mortality of the different prey groups. Thus, the crabs’ biomass and catch trends in the Visayan Sea manifest a combination of biological and behavioral responses to the changing environmental conditions, the increasing fishing effort aimed at blue swimming crabs, and a trophic response to not only reduced predation, but also reduced competition in the system as well.

Chlorophyll a and primary productivity were shown to improve the fits not only of the plankton-feeding small pelagic fish, but also of the demersal fish and invertebrate groups. This reiterates our earlier contention that, in terms of energy utilization within the Visayan Sea, there is a lack of a clear-cut delineation between the pelagic and demersal groups. Such a scenario echoes the benthic-pelagic coupling mechanism (Ricci et al., 2022) wherein benthopelagic groups act as important couplers of energy (or biomass) in the system. In the shallow Visayan Sea, zooplankton, invertebrate groups (such as squids and shrimps), and even the juvenile stanzas of the pelagic fish groups facilitate both the upward (i.e., from the benthic producers to the pelagic consumers) and the downward flow of energy (i.e., from the pelagic producers to the benthic or demersal consumers) as they feed and are similarly preyed upon by consumers across the entire water column.

The model estimated primary productivity anomalies were found to be correlated with two regional climatic indices, the Pacific Decadal Oscillation (PDO) and the Southern Oscillation Index (SOI). However, a stronger correlation was shown with the former, which is consistent with the findings in Bacalso et al. (2023). As an index of SST variations in the northern Pacific Ocean, positive PDO anomalies were shown to cause irregular monsoon wind patterns and regional SST measures in the western North Pacific (Wu, 2013; Yang et al., 2018), which have implications on the duration, frequency, and timing of the monsoon rains and typhoon landfall. As a non-upwelling system, chlorophyll concentrations in the Visayan Sea are attributed largely to land-based nutrient run-off (Cordero-Bailey et al., 2004), which affects the primary productivity levels that support the production of small pelagics in the area. However, rainfall levels that are above the seasonal averages, particularly during the rainy northeast monsoon, were found to have negative effects (Bacalso et al., 2023) due possibly to reduced salinity (Sunye and Servain, 2002) and increased eutrophication and anoxia (Holmes et al., 2021). These processes explain the chief influence of rainfall on the small pelagics in the Visayan Sea model simulations, which describe a different mechanism for seasonal primary productivity compared to wind-driven upwelling systems in the Philippines (Cordero-Bailey et al., 2004) and elsewhere (Cury and Roy, 1989; Sartimbul et al., 2010; Teixeira et al., 2016). Wind and the timing and amount of rainfall were likewise found to be crucial factors that affect the recruitment success of different fish and invertebrates in the Visayan Sea (Bacalso et al., 2023). Thus, although the contribution of the other 3 environmental parameters to the improvement of the overall model fit is lower than that of the SST, their contributions to refining and explaining the individual group fits of catch and biomass are not irrelevant. As shown in the visualized sardine catch trend, a major and a minor peak occurred in 2003 and 2014, respectively ( Figure 6A ). Similar patterns were displayed by the catch trends of mackerels ( Figure 6B ), small demersal fish ( Figure 6D ) and squids ( Figure 6E ). At the same time, Figure 2A shows that these were the same years when low SSTs coincided with high chlorophyll-a levels and moderate rainfall. Patterns of water mass movements in the Philippine archipelago correspond with environmental patterns associated with the monsoons, and were found to influence the distribution, abundance, and the species richness of zooplankton across different geographic regions (Noblezada and Campos, 2012). This finding is particularly relevant in the Visayan Sea ecosystem where zooplankton acts as a link for the flow of primary production to the pelagic and demersal consumers alike.

The only group for which no model was able to adequately explain the historical observations in catch and biomass were the scads. The adjusted vulnerability settings also imply contrasting trophic controls on the juvenile and adult stanzas. We suspect that there might be a need to further analyze the species and size composition of this group in the model. At present, it is composed of both small scad species (e.g., Selaroides leptolepis with a mean length of around 13 cm in the catch) and larger scads (e.g., Decapterus species with a mean length 26 cm in the catch). Segregating them into two groups with distinct inputs of vital rates, diets and catches may produce better results. On the other hand, the biomass of sardines appeared to be stable over time. Unlike the sardine catch trend, it showed no apparent visual correspondence with the environmental parameter trends ( Figure 2A ) nor with any of the effort trends of fishing gears targeting the sardines, such as purse seines (PS), mid-water trawls (MWT) and drift gillnets (DGN) ( Figures 2B, C ). This contrasts with the results in (Bacalso et al., 2023) where the sardine CPUE displayed a fluctuating and overall declining trend. In that study, however, only the Sardinella gibbosa species was represented, while several sardine species are represented in the present study. Over the years, the contributions of the different sardine species to the overall catch showed shifting dominance from S. gibbosa during the earlier years of assessment to S. fimbriata and S. lemuru in more recent years (Mesa, 2019). This may explain the apparent “stability” of the sardine biomass in the model. However, this may warrant further investigation as the different sardine species may have different population resilience to fishing and environmental disturbance. In other Philippine fishing grounds, the dominant sardine species captured also show earlier signs of over-exploitation, such as decreasing catch volumes with increasing fishing effort while the catch and population status of the less-exploited sardine species are sustained (Olano et al., 2009; Ramos et al., 2009). Unfortunately, local studies on the causes of sardine species population fluctuations other than fishing are lacking. Even with their temperate species counterparts, historical collapses and recovery that can be explained by alternating optimal and sub-optimal environmental conditions are explored in relation to other small pelagics, such as anchovies (Lecomte et al., 2004; Bonanno et al., 2016). This outcome thus underscores an important challenge in modeling complex and multi-species ecological systems, specifically in making decisions concerning the tradeoffs between achieving a certain level of realism and model parsimony (Coelho et al., 2019). Nevertheless, we believe that our results provide convincing evidence for the feasibility, applicability, and necessity of the ecosystem modeling approach to assessing the status of the Visayan Sea and the variable influence of fishing and environmental drivers.

The present study was able to show that fishing and its corresponding trophic effects are the major drivers of biomass and catch in the Visayan Sea over the past two decades, while the environmental drivers provide explanations for patterns that cannot be explained by fishing alone. Thus, we strongly recommend the use of ecosystem modeling in the assessment of exploited marine ecosystems to complement the analyses from typically single-species population assessments in the Philippines, in order to appropriately respond to the nationally-legislated and internationally-recommended application of the ecosystem approach to managing the country’s fisheries. As per convention, harvest controls are established for stocks relative to a set of biological reference points. Broad regulations in fisheries are also to be imposed when stocks are deemed overfished within a certain fishing ground based, for example, on the computed maximum sustainable yield (MSY). At present, none of the fisheries reference points in the country are derived from the basis of environmental effects and trophic interactions. Furthermore, none of them can be considered as an adequate indicator of fisheries at the ecosystem level. However, the recent delineation and designation of the Philippine fishing grounds as Fisheries Management Areas (FMAs) presents an opportunity to endorse the adoption of more holistic and system-level operational indices in Philippine fisheries assessments.

While the application of ecosystem modeling presents a promising solution to address this gap, it likewise warrants the enhancement of the capabilities of local fisheries analysts in this regard. Further, a systematic monitoring of relevant local environmental parameters needs to be established and tailored to complement the fisheries monitoring and assessments specific to the country. In addition to the 4 environmental parameters analyzed in this study, salinity, dissolved oxygen, and bottom temperature, which together are also considered as water quality indicators, may be important to consider in future analyses given that they are shown to also have significant effects on fish stock abundance and assemblage (Zhang et al., 2009; Duque et al., 2020; Chiarini et al., 2022).

A model that can reasonably replicate observed historical trends provides a powerful platform to forecast potential future effects of fisheries policies and management measures (Rose and Cowan, 2003; Coro et al., 2016; Capotondi et al., 2019). Of parallel relevance are the future scenarios of climatic change predicted by the Intergovernmental Panel on Climate Change (IPCC). In the Philippines, the average temperatures are expected to rise by as much as 0.9°C to 1.9°C (under RCP 4.5) and 1.2°C to 2.3°C (under RCP8.5) by the mid-21^st century (PAGASA, 2018). The same report describes increasing trends in annual and seasonal rainfall in many parts of the country, while longer seasons of lower-than-average rainfall are predicted in others. Further, declines in chlorophyll-a and primary production (by 0.5% to 11.8% of present-day mean values) are predicted following the projected increases in sea surface temperature and the shallowing of the mixed layer depth around the Philippines (Geronimo, 2018). Thus, the Visayan Sea stocks are expected to be negatively affected under these climate change projections, and the predicted economic impacts are likewise severe (Suh and Pomeroy, 2020). However, our results have demonstrated the modulating effects of trophic interactions to system disturbance over time. It would therefore be interesting to explore how the stocks in the Visayan Sea are predicted to respond to the various climate change scenarios using the ecosystem modeling approach, and moreover, be able to provide alternative insights to the future exploitation patterns and management scenarios for the Visayan Sea stocks (such as spatio-temporal effort regulations and other harvest controls).