Construction and Interpretation of Particle Size Distribution Spectra From 19 Ecopath Models of Chinese Coastal Ecosystems

To assess the changes that fisheries have imposed on the functioning of coastal marine ecosystems in China, 19 published Ecopath models were used to construct particle size distribution (PSD) spectra. The results show that high biomass of jellyfish from ranching operation impacted almost all of the ecosystems studied here. As well, an increasing impact of fisheries was demonstrated, via steeper PSD slopes, for ecosystems with models covering two or more periods. Models of nearshore areas, i.e., bays and estuaries, exhibited steeper PSD slopes than models of offshore areas. The PSD slopes were also correlated with total catch (TC), mean trophic level of catch, and the Shannon–Wiener diversity index (H′), which can be computed from Ecopath models. A multiple regression predicting the PSD slopes from year and mean trophic level of the catch explained 38.5% of the variance in the slopes. Overall, this study confirmed the status of depletion in China's marine fisheries resources. In addition, it established that including a PSD while constructing an Ecopath model, which is straightforward, will improve it, and allow more insights to be obtained from it regarding the impact of fishing on marine ecosystems.


INTRODUCTION
Marine ecosystems provide habitats to marine organisms and supply them with their nutrients (Odum and Barrett, 1971); however, the diversity of species and the multitude of interactions between them are difficult to identify and represent in models (Limburg et al., 2002;Fulton et al., 2003). Nevertheless, the Ecopath with Ecosim (EwE) modeling approach has become widely used, notably because of the relative ease with which its parameters can be estimated (Christensen and Pauly, 1992;Pauly et al., 2000;Coll et al., 2015;Colléter et al., 2015) and its ability to realistically represent the structure of food web within a given period (with Ecopath), then to simulate over time the changes of trophic interactions between species or group thereof (with Ecosim).
China has a long coastline, often divided into three large marine ecosystems (LMEs), i.e., the Yellow Sea (including a marginal sea, the Bohai Sea), the East China Sea, and the South China Sea (Wang and Aubrey, 1987;Sherman and Hempel, 2008). The structure and functioning of these ecosystems have been impacted by climate change (Jiao et al., 2015;Liang et al., 2018), ocean acidification (Liu and He, 2012), and especially overfishing (Wang et al., 2008;Liang and Pauly, 2017a,b;Liang et al., 2018;Pauly and Liang, 2019;Zhai and Pauly, 2019).
The depletion of fisheries resources in China's coastal seas was partly mitigated by the development of mariculture. Thus, in 2018, the total coastal areas devoted to marine aquaculture was above 20,000 km 2 , and production was more than 2 million tonnes (Anonymous, 2019). Indeed, the high density of mariculture enterprises and the escape of farmed animals exacerbate coastal pollution (Cao et al., 2007).
Another issue is worldwide proliferation of jellyfish (Brotz et al., 2012). While jellyfish are almost everywhere perceived as nuisance species (e.g., Lynam et al., 2005), several species are sought after in China, and some jellyfish are even farmed, or rather "ranched" (Dong et al., 2008). The increase in jellyfish in China's coastal seas in recent years are well-documented Dong et al., 2010) as is the negative impacts of jellyfish blooms on fishery and ecosystem functioning (Xian et al., 2005). For these reasons, Pauly et al. (2008) pointed out the important role of jellyfish in many ecosystems and the need to include them in food web models, which may have inspired the authors of subsequent Ecopath food web models to explicitly include jellyfish in their food webs (Lamb et al., 2019).
While the South China Sea, whose living resources are shared between numerous countries , was modeled in the early 1990s (Pauly and Christensen, 1993), the first Ecopath model of a uniquely Chinese marine ecosystem was that of Tong et al. (2000), which sought to represent the structure of the food web in the Bohai Sea in 1982-1983. Since this pioneering effort, all of China's LMEs have been modeled by Chinese authors, i.e., the Yellow Sea (Lin et al., 2013), the East China Sea , and the South China Sea (Cheung, 2007). However, many due to a perceived "lack of data, " none of these authors have used Ecosim, the dynamic routine of EwE.
Here, we use 19 extant Ecopath models of Chinese coastal ecosystems for inferences on the impact of fisheries on such ecosystems using particle size distribution spectra (Sheldon et al., 1972;Guiet et al., 2016), based on these 19 models, but after they were modified such as to generate clear linear spectra. The resulting particle size distributions (PSD) spectra were then used to describe and quantify the effect of fisheries and environmental effects.

Sources of Materials
The Ecopath models used in this study were sourced from the EcoBase (http://ecobase.ecopath.org/) and the scientific literature ( Table 1). The geographical distribution of modeled ecosystems ranged from northern to southern China (Figure 1) and involves major main types of marine ecosystems, i.e., LMEs, a marginal sea, shelves, gulfs and bays, and estuaries. In addition, the models cover the 1960s to the 2010s, which encompass a period of huge changes in China's economy, including fisheries.

Basic Principles of Ecopath Modeling
Ecopath was initially conceived by Polovina (1984) as a food web model with a set of functional groups linked by biomass flows according to the principle of mass balance, and further developed by Christensen and Pauly (1992) and . The model is structured by a master equation, which applies to all consumer groups in the models: where B i and B j are the biomass (here in tonnes wet weight by km −2 ) of prey i and predator j, respectively; P/B i is the production per biomass (tonnes km −2 year −1 ) of prey I; Q/B i is the consumption per biomass (tonnes km 2 year −1 ) of prey i; DC ij is the proportion of prey i in the diet of predator j; EE i is the ecotrophic efficiency, i.e., the fraction of i's production that is consumed within the system; Y i is the fishery catch (tonnes km −2 year −1 ); E i is the net migration rate (i.e., emigrationimmigration; in tonnes km −2 year −1 ); and BA i is the biomass accumulation rate (year −1 ).
Another key equation of Ecopath is where Q i is the consumption (tonnes km −2 year −1 ), P i is the production (tonnes km −2 year −1 ), R i is the respiration (tonnes km −2 year −1 ), and UA i is the unassimilated food of group i (tonnes km −2 year −1 ). The model assumes that the mass and energy input and output must be balanced within the time period covered by the model.

Particle Size Distribution
Particle size distribution (PSD) or size spectra within aquatic ecosystems are a powerful tool to show how much human activities and environmental-driven factors impact on aquatic ecosystem (Blanchard et al., 2017) and are usually based on sampling plankton (Hunt et al., 2015;Wallis et al., 2016) or fish (Boudreau and Dickie, 1992;Bianchi et al., 2000) using net gears that are adequate for these groups only. Irrespective of the way they are presented-which is usually as log(biomass) vs. log(body weight)-the resulting spectra frequently cover only plankton and/or or fish while regularly omitting marine mammals, i.e., they mostly cover a relatively narrow range of sizes, and thus their slopes are probably biased downward. Pauly and Christensen (2002) demonstrated that Ecopath models, which necessarily include small and large organisms with biomasses that are realistic, and mutually compatible, can be used to construct PSD, whose slopes can be compared between systems and/or with the slope of PSD estimated through field sampling. In order to be able to use the groups' biomasses of an Ecopath model in a size spectrum, two issues must be considered: (1) Many groups (e.g., large teleosts, which start as planktonic larvae) have ontogenies that span several (log) sizes classes or "bins, " and thus, their biomass must be allocated to several bins in accordance to their relative abundances.
(2) Different size/age groups (e.g., in large species) have vastly different growth rates, as do different species occurring in the same ecosystems, and thus, they remain in bins of given sizes for periods that vary among and between the species that define a group. Version 6.6 of EwE includes a routine based on Pauly and Christensen (2002), which account for these two items using, for each of the groups, the parameter of a von Bertalanffy growth function (VBGF) to express the growth of the animal therein across a range of size bins. For length, the VBGF is where L t is the mean length at age t of the species in question; L inf is the asymptotic size, i.e., the mean length attained after an infinitely long time; K is growth coefficient (here in year −1 ); and t 0 (usually negative) is the age they would have had at a size of zero if they had always grown in the manner predicted by the equation (which they have not; see e.g., Pauly, 1998). However, given that we are dealing with biomass, the VBGF here is combined with a length-weight relationship (LWR) of the form W = a·L b to obtain a version of the VBGF suitable for growth in weight, i.e., where W inf is the mean weight attained after an infinitely long time, and the other parameters are as defined previously.
Recall that, in fisheries science and in Ecopath, it is expressed by where N t is the number at time t, N t− t is the number at time t -t, and Z is the instantaneous rate of total mortality, which happens to be equal to the P/B ratio, a parameter of Ecopath models, when growth follows the VBGF (Allen, 1971). Thus, the biomass contribution B t of a species (or group, see below) in a given size bin can be computed from where t is the time the group takes to grow through a bin, or weight class, with W t as defined above; B t is scaled over all weight classes so as to sum up to the total biomass of the group. As for t, it is estimated as the time difference between the relative ages (t ′ ) computed from the inverse of the VBGF, i.e., where t ′ is the relative age corresponding to the lower or upper limit of the weight interval, and all other parameters are as defined previously. Note that t 0 , which allows absolute ages to be computed, is not required because t is computed as a difference between relative ages. Thus, to construct a PSD when a parameterized and balanced Ecopath model is available, all that is needed are L inf (or W inf ) and K values to deal with the ontogeny of each group (when L inf only is available, a LWR must be also provided, but its exponent b can be assumed = 3; Froese, 2006;Hay et al., 2020).
The growth parameters of the species reported or assumed to contribute most of the biomass in each group were used to describe the growth of each functional group. These growth parameters originated from FishBase (www.fishbase.org), SeaLifeBase (www.sealifebase.org), and the scientific literature. When several sets of growth parameters were available, those originating from China or nearby countries were selected. For species without growth parameters, similar size species from the same genus or family were taken as substitutes (see Table SI). This applied also to zooplankton, of which an increasing number of species have their growth described by the VBGF (see, e.g., Wang et al., 2017).
For Ecopath models in which jellyfish were listed as a separate group, the input biomasses and catches of jellyfish were adjusted by reducing their values to the ones they would have if they had the same water content as fish, i.e., from an average of 98% to an approximate value of 75% . These and other minor adjustments of the 19 Ecopath models (see Table 1) were done for each model separately, such as avoid artificially introducing common features.

Multivariate Analyses
Nine of the indicators output by Ecopath were used for comparisons between the 19 models at hand, i.e., the slope of PSDs, mean year of coverage, mean latitude (ML) of area covered, total system throughput (TST, in tonnes km −2 year −1 ), mean trophic level of the catch (MTL c ), mean trophic level of the ecosystem (MTL e ), total biomass (TB), total catch (TC), and Shannon-Wiener diversity index (H ′ ). MTL and H ′ values were calculated based on functional groups instead of species. All indices except mean year and mean latitude were calculated by the program EwE 6.6 .
Among these nine indicators, year, slope, MTL c , MTL e , and TC are the factors that can be related to fishery impacts on ecosystems. Four indicators, ML, TST, TB, and H ′ are environmental factors; TST is the total annual biomass flux through all direct trophic interactions, which indicate the size of ecosystems; TB and H ′ define the size and complexity of ecosystems from different perspectives; and ML deals with the impacts on ecosystems of different geographic locations (e.g., via temperature). These are assumed to be the most representative indicators of impacts on the ecosystems by both human and environmental factors.
Regression analysis was used to explore the relationships between the slopes of the PSD, and these indicators and multivariate analyses were performed to compare the structure and function of the ecosystems represented by the 19 models. Multidimensional scaling (MDS; Borg and Groenen, 2003) and H-clustering (Johnson, 1967;Bishop and Tipping, 1998;Fraley, 1998) were used to classify the 19 models based on the 9 indicators above, and the differences between the groups were tested by an analysis of similarities (ANOSIM). Similarity percentage analysis (SIMPER) was applied to identify the indicator associated with the difference between models. These analyses were conducted using Microsoft Excel and the Rstudio software.

Adjustments to the Ecopath Models
Fifteen models contained jellyfish as explicit group. They were all adjusted for water contents, i.e., there biomass was reduced (see above). In addition, their estimates of ecotrophic efficiency (EE) was increased from an average of 0.37-0.77, as suggested by a negative correlation (p < 0.01) between EE and biomass (see Table 3). This implies that their production is consumed at about twice the rate assumed by the original authors. Figure 2, pertaining to the Bohai Sea, illustrates the model adjustments that were performed to obtain the PSD slope summarized in Table 2, while the adjustments required for the other models are illustrated in Figures SI, SII. For eight of the models, the adjustments lead to the PSD exhibiting markedly  closer correlation, while the remainder of the models was not affected much. The slopes of 19 models in Chinese coastal seas are shown in Figure 3. The mean slope and standard error were −0.89 ± 0.069, and their range was −0.38 to −1.61. One of the main results is that, for all models representing 2 or more years of the same ecosystems, there was a clear steepening of the slopes with time (Figure 4), suggesting that, in recent years, the upper higher trophic levels of costal food webs have disappeared.

Multivariate Analyses
The MDS analyses and the H-clustering produced similar results, i.e., they both divided the 19 models into 5 clusters (Figure 5). One cluster consisted of northern South China Sea (1970s and 2000s) and another three models consisted of Jiaozhou Bay (1980Bay ( , 2011Bay ( , and 2015. The "stress" of MDS is an indicator of the fit of the graph to the available data. Here, the MDS result had an "excellent" fit, with a stress value of 0.01(<2.5%; Kruskal, 1964).
The ANOSIM yielded a global R of 0.70 > 0 (p < 0.001), implying that the difference between groups was significant. The  SIMPER analysis indicated that the main reason for the difference between the models was TST, i.e., the scale of the ecosystems; the second reason is TC. The slopes were significantly correlated with TC, H ′ , and TC × MTL c when all 19 models were included (see Figures 6A-C). In addition, the slopes have significant relationships with year and MTL c for the 14 models representing the same areas in different years (see Figures 6D,E). A multiple regression predicting the PSD slopes from year and mean trophic level of the catch had a multiple coefficient of variation (R 2 ) of 0.385, implying that year and mean trophic level explain 38.5% of the variance in the slopes (see Figure 6F).

DISCUSSION
The key insights emanating from this study are, in our opinion, that: (1) Constructing PSD plots from balanced Ecopath models is straightforward; (2) but jellyfish and mariculture require special treatment; (3) PSD do reflect fishing pressure on ecosystems, and thus, (4) PSD should become part of the completion and verification of Ecopath models.
We elaborate on these four points. The routine of EwE 6.6 that can be downloaded from www. ecopath.org which allows for deriving PSD from Ecopath models is, as mentioned above, very easy to use. Particularly useful is a graph that shows the biomass that each group contributes to different body weight bins and which was used for the model adjustments described above.
Moreover, contrary to the situation prevailing only a few years ago, growth parameters and LWRs are increasingly available for fish in FishBase (Froese and Pauly, 2019;Hay et al., 2020), SeaLifeBase (Palomares and Pauly, 2019), and online literature (e.g., Palomares and Pauly, 2008) such that even groups for which other growth models had been used now have the VBGF to describe their growth, a trend likely to continue.
The results obtained here suggest that one of main reasons for imbalance in Chinese coastal ecosystems is due to jellyfish, whose biomass must be adjusted-if it was no done in the original models-for their extremely high water content  and the short duration of the period during which their biomass is high-particularly in jellyfish species that are ranched (You et al., 2007;Dong et al., 2008). In these species, the high biomass period may last only 1 or 2 months, and thus, their ecosystem impact (i.e., food consumption) must be reduced accordingly; such temporal adjustments were not done by the authors of all models.
Although not used here, temporal adjustments may also have to be performed in some cases for short-lived farmed species, e.g., those on artificial reefs, deployed along coastline to simulate the natural reefs environment, to provide a replacement for lost natural habitats to marine fisheries resources (Baine, 2001). Thus, the models for artificial reefs in Laizhou Bay showed serious impacts on the trophic flows by individual farmed species, notably by a high biomass of sea cucumber and oyster that had spilled over from a nearby farm (Xu et al., 2019). Such adjustment can be done by multiplying the mean biomass of the group while it is in the system by the fraction of the year (often <<1) that it is in the system.
Another adjustment that might have to be considered in some cases (though it was not used here) is to reduce the biomass in the smallest size group (phytoplankton and microzooplankton), when its size range is smaller than the bin size of the PSD plots (such adjustments are required for reason of arithmetic and do not imply that the Ecopath models in question are in error).
As expected, different models covering the same area exhibit PSD slopes that became steeper over time, implying more fishing pressure on these areas, which is in line with what is known of fisheries along China's coastline (Ling et al., 2006;Li et al., 2015;Teh et al., 2017). The PSD for the Bohai Sea, Jiaozhou Bay, Pearl River Estuary, Beibu Gulf, and other coastal and nearshore ecosystems had steeper slope than the PSD for the wider East China Sea and the northern South China Sea, which is consistent with the result of Liang and Pauly (2017b), who demonstrated a strong fishing down effect (which induces strongly negative PSD) along the coast of the East China Sea, followed by the development of an offshore fishery, which induces another fishing down trend.
The TST index was the main reason why the northern South China Sea model is different from others, while the TC index was the main reason why the models of Jiaozhou Bay (1980Bay ( , 2011Bay ( , and 2015 formed a distinct group. Jiaozhou Bay is the most important shellfish mariculture site in the northern China, and its high TC values are mostly due to shellfish farming and fishing for escapees from farms. Therefore, it appears that mariculture has more impacts on the ecosystems than (over)fishing.
The steepness of the PSD slopes, which expresses fishing pressure, was significantly negatively correlated with the year and TC variables and was significantly positively correlated with H ′ and MTL c . This implies that the marine ecosystems that have a high biodiversity and MTL c would be in better shape than the others. Given a background of generalized overfishing, this would be reversed with higher TC, indicating that reducing the fishery catch is conducive to the balance of the marine ecosystems' structure and function. Moreover, the impacts of year variable is not direct, rather mainly associated with TC and the depletion of fishery resources, which increase through time. Thus, in China, if the high fishing pressure is maintained, or even increase, the PSD slope of coastal marine ecosystems will continue to decline. In addition, the slopes were positively correlated with MTL c , implying that an ecosystem with a less negative PSD slope will have more high-trophic level fish that can be caught, which is another definition of "fishing down" (Pauly et al., , 2001Liang and Pauly, 2019).
High level of diversity is conducive to the stability of ecosystem (Elmqvist et al., 2003), and thus, models should represent that biodiversity ought to be represented. Several of the models considered here lacked large marine mammals, sharks, and/or seabirds; while many of these larger animals do not occur along China's coast, some do, if in small numbers, and models should include them. Some PSD slopes were too steep because of the omissions.
Overall, this contribution, besides further confirming that China's coastlines are strongly exploited by fisheries, also confirmed that PSDs allow for numerous insights and inference on food web models and the ecosystems they represent. Moreover, PSDs, which are easily constructed, help identify problems with the parametrization of Ecopath models. Thus, while Ecopath-derived PSDs are not common so far (but see Palomares et al., 2009), we suggest that they should be, in China and elsewhere.

AUTHOR CONTRIBUTIONS
LZ was responsible for the data collection, formal analysis, and writing of the original draft. DP was responsible for the conceptualization, methodology, and supervision.

FUNDING
LZ research was funded by China Scholarship Council (CSC). DP research was supported by the Sea Around Us, which receives funding from the Oak Foundation, the Marisla Foundation, the Paul M. Angell Family Foundation, the David and Lucile Packard Foundation, the Minderoo Foundation, and the Bloomberg Philanthropies through RARE.