The YS is one of marginal seas of northwestern North Pacific, adjacent to the Bohai Sea in the north and the East China Sea in the south (Figure 1). The YS covers an area of about 400,000 km² with an average depth of 44 m. The YS is highly productive, contributing to more than 20% of marine fishery production in China and comprised of various species associated with different niches. Under the dual pressures of climate variability and overfishing, species composition in the total catch have changed greatly (Ma et al., 2019). In particular, the proportion of warm-water species increased significantly, while the proportion of cold-water species decreased significantly (Liang et al., 2018). In addition, the fishery species composition changed from demersal fish at high trophic levels to small pelagic fish at low trophic levels (Ma et al., 2018).

In this study, we used a series of selection criteria suggested in Otto et al. (2018a) to select candidate indicators. As fishery-independent data from research surveys are rare in the YS (e.g., Liang et al., 2020), we took advantage of the YS catch time-series data from 1960 to 2014 obtained through the Sea Around Us project¹ (Pauly and Zeller, 2016) to derive ecological indicators. Catch data were summarized to species, genus, family, order, and class. All species were classified into four functional groups, namely the large predatory, small pelagic, demersal, and invertebrate groups according to their trophic levels, biological characteristics, and previous studies (Tian et al., 2006, 2013). All the species were also categorized into three thermal groups, namely warm-water, temperate-water, and cold-water groups according to their optimal temperature (Ma et al., 2019). Based on the catch data of the corresponding groups, we calculated 22 IFCSs that represent the YS large marine ecosystem (Table 1). These 22 IFCSs include 5 indicators based on all species (termed “total-IFCSs”), 11 indicators based on functional groups (“functional-IFCSs”) and six indicators based on thermal groups (“thermal-IFCSs”). Specifically, total-IFCSs include total catch (ToC), the mean trophic level (MTL), the mean temperature of catch (MTC), Whilm species diversity (H″), and McNaughton dominance (Dom). Functional-IFCSs contain catch of large predatory groups (LPC), catch of small pelagic groups (SPC), catch of demersal groups (DeC), catch of invertebrate groups (InC), the ratio of LPC to ToC (LPC/ToC), the ratio of SPC to ToC (SPC/ToC), the ratio of the ratio of DeC to ToC (DeC/ToC), the ratio of InC to ToC (InC/ToC), the ratio of DeC to SPC (DeC/SPC), the ratio of LPC to SPC (LPC/SPC), and functional evenness based on functional groups (F-J′FD). Thermal-IFCSs consist of catch of warm-water groups (WWC), catch of cold-water groups (CWC), the ratio of WWC to ToC (WWC/ToC), the ratio of CWC to ToC (CWC/ToC), the ratio of WWC to CWC (WWC/CWC), and functional evenness based on thermal groups (T-J′FD).

Prior to the calculation of each indicator, multiple imputation with Bayesian linear regression techniques (Buuren and Groothuis-Oudshoorn, 2011) was used to fill the missing values of the catch data in the YS (1.43%). Multiple imputation was performed using the packages “mice” within the “R” statistical and programming environment (R Core Team, 2019), with each imputed value derived from five iterations and each iteration consisting of 100 imputations.

A total of eight large-scale climate indices were tested for their effects on all IFCSs, including the Pacific Decadal Oscillation (PDO), North Pacific Gyre Oscillation (NPGO), Arctic Oscillation Index (AOI), Pacific-North America Index (PNA), North Pacific Index (NPI), Southern Oscillation Index (SOI), Asian Monsoon Index (MOI), and Sea Surface Temperature in Nino 4 Region (Nino 4). Descriptions, data sources, and spatial-temporal resolutions of these climate indices are provided in Supplementary Table 1. These climatic variables are well documented and largely associate with the ecosystems and their fish communities in the North Pacific (Boeing and Duffy-Anderson, 2008; Ceballos et al., 2009; Litzow et al., 2014; Kidwell et al., 2017).

In addition to the large-scale climate indices, seven well-documented environmental indices in the YS, including sea surface temperature (SST), sea surface salinity (SSS), sea level pressure (SLP), scalar wind (SW), wind direction (WD), air temperature (AT), and precipitable water content (PWC), were also selected to represent regional environmental conditions (Cui and Zorita, 1998; Park et al., 2015; Zheng et al., 2015) and to test the responses of IFCSs. More details on these environmental indices are provided in Supplementary Table 1.

Anthropogenic pressures in the YS were approximated by fishing efforts (FE) in the following statistical analyses to understand the responses of IFCSs. Engine power data for the total number of Chinese marine fishing boats in the YS were derived from Chinese Fishery Statistics from 1960 to 2014 to represent fishing efforts (Zhao et al., 2015). Due to the lack of other accurate fishing effort indices for specific species/groups, we assumed fishing efforts were the same across all fish species.

The indicator-testing framework of Otto et al. (2018a) was used to assess the performances of IFCSs and further choose a set of complementary and robust indicators for representing the FCS in the YS. Before conducting statistical analyses, all the 22 IFCSs were standardized by zero-mean normalization. All calculations were done using the package “INDperform” (Otto et al., 2018b) within the R environment (R Core Team, 2019). The indicator-testing framework, combined with score criteria (Supplementary Table 2), is a seven-step process with the first five steps being used to validate and quantify the performances of IFCSs in their responses to pressures (Supplementary Table 1).

Specifically, Step 1 identifies indicators’ trends and links to the score criterion 1. Generalized additive models (GAMs; Hastie and Tibshirani, 1990) are applied to model the long-term changes of IFCSs in response to one pressure (explanatory variable) at a time. When there are significant temporal auto-correlation in residuals, the alternative Generalized Additive Mixed Models (GAMMs) are used. A score of 1 is given if an IFCS shows a significant trend, otherwise, score 0. Step 2 identifies the responses of IFCSs to pressures and links to sub-criterion 2.1 and the first two sub-criteria of criterion 3. Step 3 measures the robustness of an IFCS in relation to one specific pressure based on the degree of fit for the last 3 years between predictions from the model using the training dataset (i.e., excluding the last 3 years) and the test dataset (i.e., the last 3 years). The score of sub-criterion 3.3 is according to the value of the normalized root mean square error (NRMSE) on the test dataset.

In the case of non-linear responses of IFCSs to pressures, the method of finite differences (Trenkel and Rochet, 2009) is used to test the first derivative of the smoothing function, representing the instantaneous change rate of IFCSs to that of the pressures in Step 4. Confidence intervals (CI) of the first derivative are generated by carrying out a conditional bootstrap by resampling from the GAM (or GAMM) residuals. The IFCSs’ responses to changes in pressures are not considered to be statistically significant when zero is contained within the CI of the first derivative. Subsequent scoring sub-criterion 2.2 was done via counting the proportion of points inside the CI. In Step 5, the threshold GAM (TGAM; Ciannelli et al., 2004) is used to detect interactions between pressures corresponding to criterion 3.4. A zero score for sub-criterion 3.4 is given if interactions with other pressure(s) are detected. The IFCS with summed scores of sensitivity and robustness > 50% of maximum score is considered to have good performances in their responses to pressures.

In Step 6, Nightingale rose diagrams are used to visualize the overall performance of the IFCSs. Finally, in Step 7, a hierarchical group-average cluster analysis (Legendre and Legendre, 1998) is carried out based on the sub-criteria scores of sensitivity (criterion 2), robustness (criterion 3), and the other two criteria (criteria 1 and 4) to recognize redundant IFCSs. Then a final indicator suite is identified based on the following criteria (Otto et al., 2018a): (1) responsiveness to at least one pressure, (2) non-redundancy, i.e., the IFCSs is not grouped with other IFCSs at the lowest level in the cluster analysis, or having higher scores in sensitivity and robustness than their redundant counterpart in the same cluster, (3) good performance (the sum scores of sensitivity and robustness >50% of the maximum score), and (4) complement, i.e., the IFCS covers the full types of pressures together with other IFCSs.

The IFCSs in different categories tended to have differing responses to pressures (Table 2). All total-IFCSs were insensitive to climate indices while being highly affected by fishing efforts. The functional-IFCSs responded selectively to fishing efforts but not to environmental indices showing primarily low performances. By contrast, the thermal-IFCSs, especially the index T-J′FD, showed high sensitivity in their responses to all three types of pressures. All IFCSs responded significantly to at least one pressure except Dom and H″ (Table 2) with only H″ showing no trend (Supplementary Figure 1). However, none of the IFCSs showed a direct response to the regional environmental index (WD) and three large-scale climate indices (NPI, MOI, and Nino 4). Therefore, these four pressure indices were excluded from the subsequent analysis. Among the 22 IFCSs that responded to one or more pressures, half showed good performances in their responses to pressures (Table 2). In terms of sensitivity, a majority of IFCSs (18 out of 22) showed good performance with scores >50%. By contrast, none of the IFCSs scored more than 50% of the total score in terms of robustness, which resulted in weak performance for many IFCSs (Table 2). Among all IFCSs, the indicators of LPC/ToC, CWC, and WWC/ToC showed significant linear responses to all their pressures. In contrast, relevant non-linear responses to all their significant pressures were found for SPC/ToC, DeC/ToC, InC/ToC, DeC/SPC, LPC/SPC, and F-J′FD (Table 2 and Supplementary Figure 2). Another striking feature across all IFCSs was that one or more threshold interactions were found for every IFCS (Supplementary Figure 2). For example with T-J′FD, there was a threshold interaction with PNA: under lower PNA conditions (<0.97) the negative effect of fishing efforts remained constant and linear, but under higher PNA conditions (>0.97) there were no significant responses (Supplementary Figure 2).

Based on the performances of the 22 IFCSs and the need to reduce redundancies found amongst the 22 IFCSs (Figure 2), we identified two combinations (Combinations I and II) with each consisting of five indicators. Either combination can be considered as a final robust indicator suite. Combination I comprised one thermal-IFCS (T-J′FD), one functional-IFCS (LPC/ToC) and three total-IFCS (ToC, MTL, MTC). Combination II consisted of one thermal-IFCS (T-J′FD), two functional-IFCSs (LPC/ToC, SPC) and two total-IFCSs (ToC, MTL). Together, these IFCSs would allow the evaluation of impacts from large-scale climatic and regional environmental changes as well as fishing effort. The total-IFCSs are suitable for detecting the impacts of fishing efforts and climate variability with ToC and MTC showing strong responses to PDO and PNA, respectively. By contrast, MTL showed a specific response to fishing efforts. As a synthetic indicator of thermal-IFCSs, T-J′FD showed a strong and robust response to climate indices (i.e., PDO and AO), environmental indices (i.e., SSS and SW), as well as fishing efforts. For the functional-IFCSs, LPC/ToC only responded to environmental indices, but not to fishing efforts and climate indices, while SPC showed strong response to fishing efforts and two climatic indices (Table 2).

Both robust combinations of IFCSs showed that the FCS in the YS has departed substantially from their reference state (Figure 3). Although the Euclidean distance scalar of combination I suggested the YS returning to the earlier state after 2000 (Figure 3A), the state approach based on the PCA indicated clearly that the last five years of the current period are outside the reference domain in terms of PC1 dimensions (Figure 3B). There were also true for Combination II (Figures 3C,D). The first two PCs accounted for 85.5 and 88.5% of the final IFCSs variances for the Combinations I and II, respectively (Figures 3B,D).

For Combination I, step-like changes occurred in 1971/72, 1985/86, and 2003/04 for PC1 and 1976/77 and 1990/91 for PC2 (Figures 4, 5); similar results were found in Combination II (Supplementary Figures 3, 4). As all five robust IFCSs have high loadings on one or both PCs (Figure 6 and Supplementary Figure 5), this can be interpreted as a deviation from the reference period in all components of the FCS represented by the suites with values of ToC and MTC (alternatively SPC) being lower while MTL, LPC/ToC and particularly T-J′FD being higher during the reference period (Figure 6 and Supplementary Figures 1, 5).

In general, the first five steps of the indicator-testing framework have indicated that the majority of the 22 IFCSs of the YS showed good performances in terms of their sensitivity to pressures, but weak performances in terms of robustness. Regime shifts of marine fish community are frequent occurrences accompanying climate regime shifts (Tian et al., 2006, 2008; Litzow and Ciannelli, 2007; Möllmann and Diekmann, 2012), and such abrupt changes can only be explained by non-linear state changes in response to sudden changes in pressures (Petersen et al., 2008; Fu et al., 2015). However, the IFCS-pressure relationships in our study were largely linear (27 of the 51 significant pressure responses), which contribute to low scores in one of the robustness criteria i.e., Criterion 3.2. Such linear relationships may partly be due to the standardization of data as suggested by Otto et al. (2018a). In addition, all IFCSs that had significant responses to at least one pressure were found to have one or more threshold interactions, which was a primary reason for the low scores in the robustness aspect. Threshold interactions are often used to test non-additive interactions between two types of pressures, which are important for assessing population dynamics (Ciannelli et al., 2004). Threshold responses in potential interactions between pressures provide guidance for potential reference points for fisheries management (Fu et al., 2019). From a conservation perspective, a threshold may represent a risk point where the risk of species loss increases sharply (Toms and Villard, 2015). For example, overfishing cod (Gadus morhua) may change the dynamics of predation-prey relationships between cod and zooplankton, and thus damage the population dynamics of cod (Casini et al., 2009). Therefore, thresholds may offer critical insight and useful guidelines when developing ecosystem-based conservation objectives (Samhouri et al., 2010). From this point, IFCSs that show threshold interactions in their responses to pressures proved their utility for EBFM.

Although 22 IFCSs were considered from different population perspectives, only five indicators were contained in the final indicator suite which were identified based on the performance of each IFCS and the need to reduce redundancies, which were ubiquitous among the IFCSs derived from the unitary data (i.e., only catch data). The five indicators within Combination I, the thermal-IFCS T-J′FD jointly with the functional-IFCS LPC/ToC and the total-IFCSs MTL, ToC and MTC, would allow the evaluation of management measures related to all significant pressures. In particular, the introduction of comprehensive indices MTL and MTC helped improve the data dimension (going beyond catch data only) by incorporating biological characteristics of fish. The MTL and MTC indicators have been widely used in related research to reflect the trophic level and optimal temperature of species, respectively (Tsikliras and Stergiou, 2014; Hermida and Delgado, 2016). In Combination II, SPC was included in the place of MTC as both of them responded to the same pressure PNA with the same performance. Although the results of the two combinations are similar, we focused on Combination I in our result presentation and discussion as the final indicator suite to reflect long-term trajectory of the FCS in the YS as the indicator MTC incorporates biological characteristics of fish while SPC is unitary (based on catch data only). To facilitate the move towards EBFM in the YS by employing and evaluating more comprehensive ecological indicators in terms of their responses to management-related key pressures, we advocate continuous and consistent research surveys to collect pluralistic data, such as length (Large et al., 2013), abundance (Shannon et al., 2009), and age structure (Shin et al., 2010), which are largely lacking in the YS.

In previous studies, LPC (Myers and Worm, 2005) and DeC (Newman et al., 2015, 2016) were found most suitable and subsequently implemented in fisheries management. However, these indicators did not perform necessarily best for the study period (1960–2014) in the YS. Specifically, LPC had high redundancy with ToC and T-J′FD within the final indicator suit in terms of its response to fishing efforts and PDO. In addition, the relationship between DeC and its significant pressures was less robust than LPC/ToC. Such results demonstrate that not all indicators are suitable for fisheries management in all ecosystems, even if they are sensitive in their responses to relevant pressures. We would encourage that the evaluation of indicator performances is done on an ecosystem-specific basis because of the unique characterizes of an ecosystem and the diverse fishery exploitation history (Fu et al., 2019).

Furthermore, the final indicator suite would not be applicable for evaluating the impact of four pressures (i.e., WD, NPI, MOI, Nino 4) as all IFCSs were not found significant in their response to them. Therefore, more comprehensive ecological indicators are needed to evaluate the impacts of these four pressures, and the performance of candidate indicators should be validated before they can be used for assessing ecosystem dynamics and guiding fisheries management.

Among multiple pressures that impact marine ecosystems, climate change plays one of the most important roles (Barange et al., 2014; Kirkman et al., 2015; Kadin et al., 2019). For instance, climate change was found to be important factor in evaluating the responses of small pelagic fish, an important functional group in marine fish communities, to fishing pressure and decadal variations (Tian et al., 2013; Ma et al., 2018). However, climate change was not identified as important in our study. The SPC was shown to have significant response only to one climate index (i.e., PNA) with weak robustness contrary to previous studies. Moreover, some indicators such as MTC commonly used to reflect climate changes (Tsikliras and Stergiou, 2014) showed no sensitivity to any climate indices. As a whole, 13 of the 22 IFCSs tested significant in their response to one or more climate indices (Table 2). The response mechanism of fish communities to climate change is complex. Climate change may impact the dynamics of marine ecosystems via complex “Atmosphere-ocean-ecosystem” processes (Ma et al., 2019), causing increases in sea water temperature, which further alters reproduction effectiveness (Palomera et al., 2007), survival (Martins et al., 2012), distributional shifts along latitudinal and depth gradients (Perry et al., 2005; Simpson et al., 2011), and other behavioral variations (Goniea et al., 2006; Audzijonyte et al., 2020). It is generally acknowledged that the responses of marine fish communities to climate change can vary among regions (Barange et al., 2014; Guenette et al., 2014), resulting in regional specificity in the performances of IFCSs. Again, our results advocate that it is imperative to validate indicators in different ecosystems before using them in management decisions, as suggested by Östman et al. (2017).

In general, IFCSs have different sensitivity in their responses to different environment variables, and significant responses to regional environment variables mainly occurred in functional-IFCSs and thermal-IFCSs (Table 2). Regional environmental conditions affect FCS directly by limiting the distribution of organisms within their tolerance range or indirectly by altering species interactions (Ochoa-Hueso, 2016). For example, Ocean warming has modified marine biology community structure by increasing biodiversity and miniaturizing the plankton ecosystem, which negatively affects demersal Atlantic cod (Beaugrand et al., 2010). Younger copepodites have strong sensitivity to midwater temperature in spring and summer, while there is a strong positive correlation between females and older copepodites and deep water salinity from spring till autumn (Otto et al., 2020). Moreover, potential threshold interactions among pressures can also alter the shape of the relationships between responses and pressures (Pang et al., 2018). Therefore, considering species compositions in specific food webs is important when employing ecosystem indicators in the evaluation of pressures and management measures.

Marine ecosystems are affected by multiple human pressures (Halpern et al., 2008). Early fisheries studies around the world have provided sufficient proof that several fish groups have declined markedly in population sizes even resulting in local extinction of some vulnerable species, both as target and by-catch species (Stevens et al., 2000). Continuous overfishing also resulted in miniaturization of important commercial species and reconstruction of total fish communities (Conover et al., 2005). Consequently, fishery-based ecosystem indicators are commonly sensitive to fishing efforts (Gascuel et al., 2016). Although we only used engine power data as proxies for the total number of Chinese marine fishing boats in the YS, our results showed that 15 of 22 IFCSs were sensitive to fishing efforts. In particular, the indicators MTL, LPC and T-J′FD were found to have highest sensitivity and robustness compared with other indicators within the final indicator suite. In contrast, the indicators DeC, LPC/ToC, DeC/ToC, DeC/SPC, WWC/ToC showed no significant responses to fishing efforts. The responses of indicators to fishing can be altered by species interactions. For example, small pelagic forage fish often increases apparently accompanying the decline of large predatory fish due to overfishing, which in turn results in the decrease of the biomass of zooplankton and phytoplankton (Frank et al., 2005). Therefore, Using fishing mortality of specific species is optimal when assessing the influences of fishing. In future studies, we intend to introduce more accurate mortality rate of specific species as proxies for anthropogenic pressures.

Combining the scores of IFCSs with the hierarchical group-average cluster analysis, we identified five IFCSs forming a final indicator suite as a proxy for the FCS in the YS. Although only 5 of 22 IFCSs were retained in the final indicator suite and the basic data were simple, our analysis sheds light on the dynamics of the FCS in the YS, proving sufficiency of the five IFCSs in the final indicator suite to represent synthetical information on fish community dynamics. Regime shifts detected in the mid-1970s and the early 1990s (Figure 4) coincided with the climate regime shifts based on the integrated trend analysis using many climate indices, while detected shifts in the final indicator suite in the mid-1980s coincided with the regime shifts of regional environment variables (Ma et al., 2019). The shift detected in the early 1970s was not accompanied by any types of pressures, which may result from the fact that no indicators tested in this study were sensitive to the four pressures (i.e., WD, NPI, MOI, Nino 4), indicating that the final indicator suite was unsuitable for evaluating the impacts of these four pressures.

Both state space approaches indicated that the status of the FCS of the YS deviated obviously from the reference year, despite some level of returning to the earlier state after 2000 (Figure 3A). There is ample evidence that the dynamics of the FCS in the YS synchronized with other regions in the world based on more relevant indices (Perry et al., 2005; Atkinson et al., 2011; Kirkman et al., 2015). The long-term trajectory of the final indicators suite in the study period derived from Euclidean distances analysis was highly consistent with fishing effort in the YS, suggesting that the final indicators suite could reflect fishing efforts appropriately. The PCA convex hull analysis for the final indicator suite showed that SPC contributed most to the observed changes in the IFCSs state space, reflecting the declines of large predatory fish biomass in global ocean (Ransom and Worm, 2003) and also reflecting the increases of small pelagic fish because of reduced predation mortality by large predatory fish (Cury et al., 2000). Moreover, the fluctuation of small pelagic fishes varied with species in the YS (Ma et al., 2018), implying that the increases of small pelagic fish were not only due to the reduction of predators. Therefore, our results also indirectly indicated that environmental conditions may have changed compared to the reference period as small pelagic fish are highly dependent on external environmental conditions (Merino et al., 2010). Such dynamics are also reflected in the decline of the mean trophic level (Zhang et al., 2007), supported by our result concerning MTL. In combination I, MTC plays the most important role, reflecting that global warming is more friendly to warm-water species than cold-water species (Almodóvar et al., 2012; O’Gorman et al., 2016). Global warming caused water temperature range to gradually fluctuate away from the niche of cold-water species, which is adverse to cold-water species but advantageous to warm-water species. Such results were also reflected in another important indicator, i.e., T-J′FD, contributing to the observed changes. Its long-term changes indicated that the resources distribution of three thermal groups trends more evenness with CWC/ToC gradually declining while WWC/ToC generally increasing within the observation period (Supplementary Figure 1). ToC ranked second in its contribution to changes in the IFCS state space, reflecting the rapid increasing of total catch (FAO, 2016).

As the YS ecosystem is continually facing impacts from overfishing and climate change, resulting in serious declines of fishery resources and miniaturization, it is imperative to move towards EBFM. The great change of the FCS in the YS has been proven by our final indicator suite, which can subsequently be used as an effective tool in ecosystem monitoring to reach management targets (Levin et al., 2009). In particular, the increase of SPC has been identified as one of the most outstanding changes of the FCS in the YS. Therefore, we advocate that the management objective of increasing size of fish should be include in the EBFM in the YS, which can potentially be achieved by enlarging mesh size.